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ABSTRACT 


The  deep  ocean  has  a  refractive  waveguide  velocity  structure  created  by 
variations  in  temperature,  salinity  and  pressure.  Because  the  sound  speed  changes 
slowly  with  range  and  average  depth-dependent  profiles  are  known  from 
measurements,  an  approximate  solution  to  the  acoustic  field  in  the  waveguide  can 
be  formed  using  perturbation  techniques.  We  then  seek  to  formulate  a  direct 
inversion  algorithm  based  on  the  perturbative  solutions.  Under  certain  restrictions, 
the  unknown  velocity  perturbations  about  a  background  profile  and  the  scattered 
field  data  can  be  written  as  a  Fourier  transform  pair. 

In  this  thesis,  we  investigate  two  well-known  perturbative  solutions;  the  Born 
amplitude  expansion  and  the  Rytov  phase  expansion.  The  approximations, 
although  closely  related,  behave  quite  differently  depending  on  the  size  of  the 
perturbation,  the  distance  traveled  in  the  perturbed  media  and  the  local  field 
gradient.  Studying  the  behavior  of  the  approximations  in  the  forward  problem 
gives  an  indication  of  the  size  and  type  of  perturbations  recoverable  in  the  inverse 
problem. 

We  begin  by  investigating  the  perturbative  solution  behavior  for  simple 
velocity  structures.  In  a  constant  velocity  or  a  layered  waveguide,  the  solutions 
and  the  first  order  errors  can  be  derived  explicitly.  Although  not  representative  of 
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a  true  ocean,  the  behavior  of  the  Born  and  Rytov  solutions  for  these  s;mple 
velocity  structures  permits  us  to  qualitatively  predict  their  behavior  in  an 
inhomogeneous  ocean  environment.  Numerical  examples  using  depth-dependent 
background  profiles  confirm  that  the  range  of  validity  of  the  distorted-wave 
approximations  is  consistent  with  that  found  in  the  constant  background  test  cases. 

The  inverse  technique  consists  of  relating  the  scattered  data  and  the 
perturbation  as  a  Fourier  transform  pair.  In  order  to  derive  the  inverse  relation, 
the  phase  of  the  transform  must  be  monotonic.  Because  the  energy  propagates  via 
refraction  in  the  ocean,  this  criteria  can  not  be  met  globally.  In  fact,  the  regions  in 
the  ocean  waveguide  in  which  the  phase  is  monotonic  are  shadow  zones;  regions 
where  the  waveguide  is  not  well-illuminated.  In  a  insonified  waveguide,  however, 
the  phase  is  not  monotonic  and  the  inversion  algorithm  is  invalid. 

Although  inverse  algorithms  using  Fourier  techniques  were  found  to  be 
inappropriate  for  perturbation  reconstruction  in  the  ocean,  the  range  of  validity  of 
the  Born  and  Rytov  propagation  models  in  a  depth-dependent  environment  was 
established.  This  has  not  been  previously  done,  and,  as  a  result,  the  knowledge 
obtained  can  be  used  to  investigate  other  inverse  formulations  based  on  these 


perturbative  solutions. 
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INTRODUCTION 


In  the  last  decade,  the  inverse  problem  in  the  ocean  has  received  considerable 
attention  from  the  ocean  acoustics  and  oceanographic  scientific  communities.  In 
1979,  Munk  and  Wunsch  proposed  a  inverse  method  that  is  now  known  as  ocean 
acoustic  tomography.  This  technique  consists  of  measuring  ray  travel  times 
between  source-receiver  pairs  and  inverting  these  to  determine  perturbations  to  an 
assumed  background  profile.  Many  experiments  have  been  conducted  proving  the 
validity  and  usefulness  of  the  technique  as  well  as  its  limitations  (for  some  of  the 
earliest  results  see  Spiesberger,  Spindel  and  Metzger  1980;  Brown  et  al.  1980;  Munk 
and  Wunsch  1982;  Spiesberger  1983;  Munk  and  Wunsch  1983).  Because  the  travel 
time  equations  axe  linearized,  the  background  or  guess  profile  must  be  known  to 
within  ~2  m/s  of  the  true  profile.  The  inversion  methodology  is  in  the  class  known 
as  indirect  inversion  techniques;  a  background  profile  is  chosen,  the  forward  model 
is  run  repeatedly  and  the  results  compared  to  the  experimental  data.  In  other 
words,  the  technique  is  based  on  iteration  of  a  forward  model. 

In  this  thesis,  we  attempt  a  different  approach  and  formulate  a  direct 
inversion  algorithm.  Using  first  order  perturbation  expansions,  the  scattered  field 
data  as  a  function  of  frequency  and  the  profile  perturbation  as  a  function  of  the 
phase  can  be  directly  related  as  a  Fourier  transform  pair.  However,  because  the 
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inversion  is  based  on  the  forward  model  we  first  must  investigate  the  range  of 
validity  of  the  perturbative  solutions  in  the  ocean  waveguide.  The  theoretical 
analysis  of  the  forward  model  is  done  in  chapter  2  and  numerical  examples  are 
presented  in  chapters  3  and  4. 

We  study  two  perturbative  solutions;  the  Bom  amplitude  expansion  (Born 
1926)  and  the  Rytov  phase  expansion  (Rytov  1937;  Chernov  1960).  The 
background  field,  on  which  the  perturbative  solutions  axe  based,  is  chosen  to  be  the 
proper  mode  solution;  the  discrete  part  of  the  normal  mode  solution.  The  normal 
mode  representation,  which  contains  both  the  discrete  and  continuous  spectrum,  is 
exact  under  certain  restrictions  (see  chapter  1).  In  the  far  field,  the  continuous 
spectrum  can  be  neglected  and  the  proper  or  propagating  modes  are  the  main 
contribution  to  the  solution. 

In  order  to  obtain  simple  error  estimates  for  a  waveguide  geometry,  we  first 
investigate  the  perturbative  solutions  in  a  constant  background  environment. 
When  the  background  is  chosen  to  be  a  constant,  the  solutions  are  known  simply 
as  the  Born  and  the  Rytov  approximations. 

Although  the  maximum  deviation  in  the  sound  speed  in  the  deep  ocean  is  only 
about  10%  over  the  entire  waveguide,  the  background  velocity  structure  in  the 
ocean  cannot  be  assumed  constant.  Therefore,  we  must  choose  a  depth-dependent 
background  sound  speed  structure  to  obtain  valid  zeroth  order  solutions  in  the 
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refractive  ocean  waveguide.  In  this  case,  the  perturbative  approximations  axe 
called  the  distorted-wave  Born  (DWB)  and  the  distorted-wave  Rytov  (DWIt)  to 
distinguish  the  solutions  from  those  based  on  a  homogeneous  background. 

Once  the  forward  model  has  been  investigated,  we  have  an  estimate  of  the 
magnitude  and  extent  of  perturbations  that  can  be  recovered  in  an  inversion.  In 
chapter  5,  we  determine  the  applicability  of  a  Fourier  inverse  technique  by 
studying  the  attributes  of  the  transform  kernel  as  a  function  of  depth  and 
frequency.  Although  no  inversions  were  obtained  due  to  phase  averaging,  band- 
limiting  and  stationary  phase  points  in  depth,  the  analysis  leads  to  a  better 
understanding  of  the  type  of  direct  inversion  algorithms  that  could  be  used  to 
reconstruct  depth-dependent  perturbations. 

The  major  new  work  presented  in  this  thesis  is  the  characterization  of  the 
direct  propagated  acoustic  field  using  approximation  methods  and  the  formulation 
of  inversion  representations  using  these  approximations. 
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1.  MODELING  THE  SOUND  FIELD  IN  THE  OCEAN  WAVEGUIDE 

1.1  Introduction 

The  underwater  sound  channel  or  SOFAR  (Sound  Fixing  And  Ranging) 
channel,  discovered  in  the  early  1940’s  by  M.  Ewing  and  J.L.  Worzel  (1948),  is  the 
most  characteristic  feature  of  the  deep  ocean  at  moderate  and  equatorial  latitudes. 
This  acoustic  channel  is  formed  by  temperature,  pressure  and  salinity  variations  in 
the  ocean.  The  temperature  and  salinity  decrease  below  the  ocean  surface  causing 
a  corresponding  decrease  in  the  sound  speed.  At  approximately  1  kilometer  (km) 
below  the  surface,  a  leveling  off  of  the  temperature  decrease  and  steadily  increasing 
pressure  causes  a  nearly  linear  increase  in  the  sound  speed  as  the  depth  is  further 
increased.  Thus,  an  acoustic  waveguide  is  created  in  which  sound  can  propagate, 
through  refraction,  for  hundreds  of  kilometers. 

Although  the  ocean  has  both  regular  and  random  inhomogeneities,  this  thesis 
will  only  deal  with  regular  or  deterministic  effects;  in  particular,  propagation  in  the 
underwater  sound  channel.  The  sound  field  is  modeled  using  the  acoustic  wave 
equation  assuming  a  pressure  release  surface  and  either  a  rigid  bottom  or  an 
infinite  half  space  below  the  waveguide  bottom.  We  will  use  a  point  cylindrical 
monochromatic  source  representation  and  assume  the  medium  is  cylindrically 
symmetric.  In  general,  the  sound  speed  is  a  function  of  both  depth  and  range. 
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In  the  next  section,  we  briefly  discuss  the  attributes  of  several  numerical 
propagation  models  in  use  today  (for  a  summary  of  numerical  propagation  loss 
models,  see  DiNapoli  and  Deavenpor*  1979);  in  particular,  the  choice  of  a 
background  model  on  which  the  perturbative  solutions  are  based.  The  last  section 
of  this  chapter  illustrates  the  variation  of  the  sound  field  as  a  function  of  source 
depth  and  frequency  for  several  sound  speed  profiles.  These  plots  are  generated 
from  the  proper  or  discrete  mode  solution;  the  background  model  used  throughout 
this  thesis. 


1.2  Numerical  Models 

A  variety  of  propagation  codes  have  been  developed  to  model  the  sound  field 
in  the  ocean.  Each  is  valid  within  a  certain  range  of  chosen  medium  parameters. 
Including  even  a  condensed  discussion  of  all  the  computer  codes  currently  in  use  is 
impossible.  Instead,  we  first  briefly  describe  two  propagation  loss  models;  the 
classical  ray  method  and  the  parabolic  approximation  to  the  Helmholtz  equation. 
Next,  we  focus  on  a  third,  the  discrete  mode  solution,  which  is  the  background 
model  used  in  this  thesis.  The  rationale  for  basing  the  perturbative  solutions  on 
this  particular  model  will  be  made  clear  in  the  remainder  of  this  section. 

The  classical  ray  equations  and  the  parabolic  approximation  are  derived  based 


on  assumptions  about  the  propagation  environment.  The  ray  equations  are  formed 


from  a  high  frequency  assumption  and  by  neglecting  diffraction  effects  (Born  and 
Wolf  1965).  The  parabolic  equation,  originated  by  Leontovich  and  Fock  (1946) 
and  introduced  to  the  ocean  acoustics  community  by  Hardin  and  Tappert  (1973), 
is  based  on  a  narrow  angle  approximation.  In  both  of  these  cases,  approximate 
equations  are  derived  from  the  Helmholtz  equation  and  it  is  from  these  equations 
that  the  solutions  are  found. 

The  ray  representation  has  a  intuitive  geometric  appeal.  Rays  are  lines 
orthogonal  to  the  wavefronts  and  so  point  in  the  direction  of  propagation.  The 
rays  can  be  found  from  the  eikonal  equation  independent  of  the  amplitude.  The 
amplitude  can  be  calculated  by  finding  the  intensity  along  each  ray  assuming  no 
coupling  between  rays.  However,  the  theory  is  invalid  in  shadow  zones  (regions  in 
which  trapped  rays  do  not  penetrate)  and  near  caustics  (envelopes  of  the  family  of 
rays)  because  the  diffraction  effects  are  neglected  in  deriving  the  equations.  As  the 
range  increases,  the  caustic  regions  become  wider,  further  limiting  the  applicability 
of  ray  theory.  In  addition,  the  eikonal  and  transport  equations  are  derived  under 
the  high  frequency  assumption;  the  length  scales  in  the  problem  must  be  much 
greater  than  the  wavelength  of  sound.  By  studying  the  ray-mode  connection, 
corrections  have  been  made  to  increase  the  accuracy  of  the  ray  models  (Pedersen 
and  Gorden  1972)  or  alternately,  a  hybrid  ray-mode  propagation  model  can  be 
implemented  (Kamel  and  Felsen  1982).  Nonetheless,  in  order  to  avoid  the 
limitations  of  ray  theory,  we  turn  to  full-wave  representations  of  the  field. 
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The  parabolic  approximation  to  the  Helmholtz  equation  is  based  on  narrow 
angle  assumption.  The  limited  angular  aperture,  characteristic  of  long  range  sound 
propagation  in  the  ocean,  makes  a  parabolic  approximation  ideal  for  the  deep  ocean 
environment.  Another  advantage  of  the  parabolic  equation  solution  is  that  it  is 
numerically  solved  using  the  split-step  algorithm,  which  permits  the  solution  to  be 
marched  in  range.  However,  the  parabolic  is  exact  only  for  horizontal  rays;  any 
angle  off  the  horizontal  results  in  a  distortion  of  the  modal  phase  which  increases 
with  range  (Brekhovskikh  and  Lysanov  1982).  As  with  the  ray  equations, 
corrections  have  been  made  to  the  parabolic  which  reduce  the  phase  error 
substantially  (DeSanto,  Perkins  and  Baer  1978),  but  nonetheless,  the  solution  is 
derived  from  a  parabolic  equation  instead  of  the  elliptical  Helmholtz  equation. 

The  normal  mode  representation  of  the  field  is  an  exact  solution  of  the 
Helmholtz  equation  if  the  sound  speed  depends  only  on  depth  and  the  boundaries 
of  the  waveguide  are  planar;  parallel  to  the  range  coordinate.  In  the  near  field  of 
the  source,  the  normal  mode  solution  includes  both  the  discrete  and  continuous 
spectrum.  Away  from  the  source,  however,  the  continuous  spectrum  is  negligible 
and  only  the  propagating  modes  contribute  significantly  to  the  field  solution 
(Ewing,  Jardetzky  and  Press  1957);  accordingly,  only  these  modes  are  computed 
numerically.  As  contrasted  with  the  ray  and  parabolic  solutions,  the  exact  solution 
to  the  Helmholtz  equation  is  found  and  then  a  portion  of  the  solution  discarded  in 
the  far  field.  The  normal  mode  solution  is  accurate  for  all  frequencies  and  the 
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proper  mode  solution  becomes  increasingly  accurate  the  farther  we  are  from  the 
source  (Tolstoy  and  Clay  1985).  Because  we  are  interested  in  long  range  low 
frequency  sound  propagation,  the  proper  mode  solution  was  chosen  to  be  the 
background  model  from  which  the  distorted-wave  Born  and  Rytov  solutions  are 
constructed. 

1.3  Examples  of  Waveguide  Propagation 

In  this  section,  several  examples  of  waveguide  propagation  are  illustrated  with 
the  use  of  color  pixel  maps.  The  field  strength  in  the  waveguide  is  represented  by 
the  transmission  loss  in  decibels  (dB).  The  transmission  loss  is  defined  by 
— 20log10  |  P/P  ref\ ,  where  P  is  the  pressure  field  and  PTtj  is  the  field  strength  at  1 
meter  (m)  from  the  source. 

In  order  to  illustrate  the  focusing  effect  that  the  SOFAR  channel  has  on 
sound  propagation,  we  compare  the  propagation  loss  in  a  constant  velocity 
waveguide  to  that  of  a  bilinear  profile  with  a  sound  speed  minimum  at  1  km  depth. 
Effects  of  source  placement  in  depth  and  changes  in  the  source  frequency  are  also 
compared.  Before  illustrating  the  acoustic  field  behavior  for  different  waveguide 
geometries,  we  first  briefly  discuss  the  discrete  modal  solution. 

Recall  that  we  are  solving  the  Helmholtz  equation.  If  the  sound  speed  in  the 
waveguide  is  constant  or  at  most  depends  on  depth,  the  homogeneous  form  of  the 
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equation  is  separable.  The  range  solution  is  an  outgoing  wave  Hankel  function. 
The  vertical  wave  function,  tp,  satisfies  the  equation 


with  boundary  conditions  at  the  surface  and  bottom  of  the  ocean.  The  index  of 
refraction,  n(z),  is  given  by  c0/c(z)  and  the  reference  wavenumber  is 
k0  =2 irf/c  0;  c0  is  the  reference  sound  speed,  usually  chosen  as  the  minimum  value 
in  the  waveguide,  c  (2)  is  the  depth-dependent  sound  speed  and  f  is  the  source 
frequency.  Because  we  are  solving  a  boundary  value  problem,  there  are  non-trivial 
solutions,  V’n  only  for  a  discrete  set  of  values  of  the  parameter,  £„.  Therefore,  the 
tpn  are  the  eigenfunctions  and  the  £n  are  the  eigenvalues  of  the  problem.  As  will  be 
illustrated,  the  number  of  eigenvalues  contributing  to  the  solution  depends  on  the 
boundary  conditions,  the  source  frequency  and  the  sound  speed  structure  in  the 
waveguide. 

Figures  1  and  2  illustrate  the  transmission  loss  in  a  constant  velocity 
waveguide  for  source  frequencies  of  25  and  100  Hertz  (Hz),  respectively.  (The 
maximum  allowable  phase  velocity  is  chosen  to  be  1600  m/s.)  The  source  is  placed 
at  a  depth  of  500  m.  We  have  a  pressure  release  surface  (P  =  0  at  the  ocean 
surface)  and  a  hard  bottom  (  P2  =  0,  where  the  subscript  denotes  the  depth 
derivative). 


Figure  1 . 


Transmission  loss  in  a  constant  velocity  waveguide.  The  source 
frequency  is  25  Hz. 
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Figure  2. 


Transmission  loss  in  a  constant  velocity  waveguide.  The  source 
frequency  is  100  Hz. 
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Because  the  velocity  is  constant  over  the  full  waveguide,  there  is  no  refraction; 
the  energy  propagates  through  surface  and  bottom  reflections.  We  have  total 
internal  reflection  due  to  the  bottom  boundary  condition;  in  general,  the  angular 
aperture  extends  from  0  to  ±90°  measured  from  the  horizontal  at  the  source 
position.  As  is  evident  from  the  isolated  bright  spots,  the  energy  focuses  only 
locally. 

The  different  structure  in  Figures  1  and  2  is  due  solely  to  the  change  in  source 
frequency.  At  25  Hz,  there  are  50  modes  contributing  to  the  field,  while  at  100  Hz, 
we  have  4  times  as  many.  A  higher  source  frequency  does  not  imply  higher  angles 
but  instead,  results  in  a  larger  number  of  modes  in  the  angular  aperture.  This  is 
easily  seen  from  the  relationship  le^n]  =  £n2  ~  A„2.  I11  a  constant  velocity 

medium,  the  index  of  refraction,  ne,  is  1,  so  that  horizontal  wavenumber 
eigenvalues  are  given  by  £„2  =  Jtj,  ~  A„2.  The  A„  vertical  wavenumbers  are  given 
explicitly  by  A„  =  (2n  —  l)n/2h  where  h  is  the  width  of  the  waveguide  and  n  = 
1,2,...  If  £mm  =  0,  we  have  only  vertical  propagation  (the  takeoff  angle  is  ±90  ° ) 
and  k0  =  ( 2N-\)n/2h  so  N  =  number  of  propagating  modes  =  k0h/ir  +  1/2. 
Determining  N  for  source  frequencies  of  25  and  100  Hz  using  this  equation,  results 
in  the  above  values  of  50  and  200  modes,  respectively. 

We  now  look  at  sound  propagation  in  a  refractive  waveguide.  In  this  case,  the 
waveguide  lies  over  an  infinite  half-space  having  a  sound  speed  that  matches  the 


sound  speed  at  the  bottom  of  the  waveguide.  Therefore,  there  are  no  bottom 
reflections;  only  surface  reflections  and  refraction.  Illustrated  in  Figure  3  is  the 
sound  speed  profile  used  to  generate  the  plots  in  Figures  4,5,  and  6.  Figures  4  and 
5  are  the  transmission  loss  plots  for  source  frequencies  of  25  and  100  Hz, 
respectively.  The  source  is  placed  on  the  channel  axis  at  1  km.  The  number  of 
propagating  modes  is  found  to  be  16  (Tolstoy  and  Clay  1987)  in  Figure  4  and  64  in 
Figure  5.  In  these  examples,  the  number  of  modes  depends  on  the  change  in  sound 
speed  above  and  below  the  channel  axis,  the  distance  of  the  channel  axis  from  the 
surface  and  bottom  of  the  ocean  and,  of  course,  the  source  frequency. 

In  contrast  with  the  constant  velocity  waveguide,  the  refraction  and 


SOUND  SPEED  vs.  DEPTH 

1470  1505  1540 


Sound  speed  profile  used  in  generating  the  results  in  Figures  4,  5  and 
6.  We  assume  an  infinite  half-space  below  the  waveguide  in  which 
the  sound  speed  is  constant  and  matches  the  sound  speed  at  the 
bottom  of  the  waveguide. 


Figure  3. 


Figure  4. 


Transmission  loss  in  a  refractive  medium  generated  using  the  profile 
in  Figure  3.  The  source  is  placed  on  the  channel  axis  at  1  km  and 
the  source  frequency  is  25  Hz. 
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TRANSMISSION  LOSS  IN  THE  OCEAN 


Figure  5. 


Transmission  loss  in  a  refractive  medium  generated  using  the  profile 
in  Figure  3.  The  source  is  placed  on  the  channel  axis  at  1  km  and 
the  source  frequency  is  100  Hz. 
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Figure  6. 


Transmission  loss  in  a  refractive  medium  generated  using  the  profile 
in  Figure  3.  The  source  is  placed  at  a  depth  of  3  km  and  the  source 
frequency  is  100  Hz. 
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channeling  of  the  energy  is  clearly  seen.  Because  the  source  is  placed  on  the 
channel  axis,  a  large  portion  of  the  energy  is  propagating  at  nearly  horizontal 
angles  with  very  little  loss.  It  is  also  interesting  to  note  the  distinct  ray-like 
structure  of  Figure  5  as  compared  to  Figure  4.  As  is  well-known,  a  ray  is  the  locus 
of  neighboring  modes  propagating  in  phase  (Officer  1958;  Tolstoy  and  Clay  1966), 
and  the  higher  the  frequency,  the  greater  the  number  of  modes.  Therefore,  we 
obtain  a  ray-like  picture  of  the  acoustic  field. 

Figure  6  illustrates  the  effect  of  moving  the  source  below  the  channel  axis;  the 
source  is  now  located  at  a  depth  of  3  km.  Although  the  waveguide  is  still 
refractive,  the  source  placement  prohibits  the  trapping  of  energy  in  the  vicinity  of 
the  channel  axis.  The  beam  has  narrowed  and  the  width  of  the  shadow  zones  has 
increased.  Instead  of  focusing  over  the  entire  range  of  propagation  as  in  Figure  5, 
the  energy  is  locally  focused  at  the  source  depth,  at  65  km  in  range. 

As  illustrated  in  the  previous  three  figures,  the  deep  ocean  waveguide  is  an 
inhomogeneous  environment  in  which  sound  propagates  mainly  through  refraction. 
This  is  contrasted  with  the  field  structure  in  a  constant  velocity  waveguide  in 
which  the  energy  propagates  via  reflection  from  the  waveguide  boundaries.  In  a 
waveguide  with  a  refractive  velocity  profile,  the  field  structure  depends  on  the 
source  frequency,  the  sound  speed  gradient  and  the  source  placement  with  respect 
to  the  channel  axis.  If  the  source  is  placed  in  the  vicinity  of  the  channel  axis  the 
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waveguide  will  be  well-illuminated.  As  we  move  the  source  away  from  the  axis,  the 
energy  is  progressively  restricted  to  a  narrower  angular  aperture  and  as  a  result, 
the  width  of  the  shadow  zones  increases. 

In  subsequent  chapters,  we  will  investigate  the  behavior  of  two  perturbative 
solutions  to  the  total  acoustic  field  in  which  the  background  solution  is  chosen  to 
be  the  discrete  part  of  the  normal  mode  representation.  Although  the  solutions 
and  their  errors  are  initially  analyzed  in  several  constant  velocity  waveguides,  we 
stress  that  a  constant  background  solution  is  not  a  valid  zeroth  order  solution  in  a 
refractive  ocean  waveguide.  The  analysis  is  done  simply  to  better  understand  and 
predict  the  validity  of  the  solutions  in  the  inhomogeneous  ocean  environment. 
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2.  THE  FORWARD  MODEL 

2.1  Introduction 

The  pressure  field,  P,  in  the  ocean,  can  be  shown  to  satisfy  the  inhomogeneous 
Helmholtz  equation.  In  general,  closed  form  solutions  cannot  be  found  due  to 
complicated  velocity  structures.  Instead,  approximate  solutions  are  constructed 
based  on  our  knowledge  of  the  characteristics  of  the  medium  in  which  the  wavefield 
is  propagating. 

One  method  used  in  analyzing  differential  equations  is  perturbation  theory. 
Perturbation  theory  permits  the  formulation  of  approximate  solutions  if  l)  a  small 
parameter,  e,  can  be  identified  in  the  differential  equation,  2)  a  zeroth  order  (c  =  0) 
solution  can  be  found  and  3)  the  true  solution  is  closely  approximated  by  the  first 
several  terms  in  the  perturbation  series  (Bender  and  Orszag  1978).  Therefore,  in 
solving  the  Helmholtz  equation  for  an  arbitrary  velocity  model  using  perturbation 
techniques,  we  must  be  able  to  construct  a  background  solution  which  is  close  to 
the  true  solution  or  equivalently,  the  guess  profile  must  be  close  to  the  true  profile. 
It  is  important  to  note  that  although  the  field  is  expanded  in  a  perturbation  series 
in  e,  this  is  for  convenience  only.  The  actual  small  parameter  is  later  identified  as 
the  difference  in  the  true  and  background  refractive  indices  and  e  wi'l  be  set  equal 


to  1. 
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Because  the  ambient  range  variability  in  the  ocean  is  at  most  a  few  percent  of 
the  depth  variability,  the  background  profile  is  chosen  to  be  strictly  depth- 
dependent.  The  pressure  field  or  zeroth  order  solution  in  this  environment  can  be 
written  as  a  modal  sum  if  the  receiver  is  in  the  far  field  of  the  source.  The 
background  solution  in  this  thesis  is  computed  using  the  Kraken  normal  mode  code 
(Porter  1985)  and  we  perturb  around  this  solution  to  form  an  approximation  to  the 
total  field.  This  is  a  major  difference  in  approach  from  that  generally  used  in  other 
disciplines  such  as  atmospheric  wave  propagation  and  seismic  exploration.  A 
constant  background  assumption  will  not  give  the  correct  zeroth  order  physics  in 
the  ocean. 

Two  perturbation  expansions  are  investigated  in  this  thesis.  The  first,  the 
distorted-wave  Born  (DWB)  solution,  is  formed  by  expanding  the  amplitude  of  the 
field,  and  the  second,  the  distorted-wave  Rytov  (DWR)  by  e'.panding  the  complex 
phase.  Although  the  DWB  and  DWR  solutions  are  derived  for  arbitrary  depth 
and  range-dependent  profiles,  only  constant  velocity  and  depth-dependent 
examples  are  investigated  numerically. 

2.2  The  Background  Solution 

The  acoustic  field  in  the  ocean  satisfies  the  inhomogeneous  Helmholtz 


equation 
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V2  +  k?0n2{r,z) 


P{ko,r,r0,z,Zo)  = 


S(r-ro)Hz~Zo) 


2n  r 


(1) 


in  cylindrical  coordinates.  P(k0,r,r0,z, z0)  is  the  total  acoustic  pressure  field  that 
we  wish  to  approximate,  V2  is  the  cylindrical  Laplacian,  n(r,z)  =  c0/c(r,z)  is  the 
normalized  sound  speed  in  a  depth  and  range  dependent  environment,  and 
k0  —  iof  c  0  is  the  reference  wavenumber.  The  source  is  located  at  a  range  r  —  r0 
and  a  depth  z  —  z0.  We  consider  the  general  case  of  r0  =*=  0  in  order  to  construct 
the  scattered  field. 


The  choice  of  the  reference  sound  speed,  c„,  is  arbitrary.  Because  the 
Helmholtz  equation  is  independent  of  e0,  so  are  its  solutions.  From  this  it  follows 
that  any  approximate  solution  must  depend  only  weakly  on  the  choice  of  the 
reference  sound  speed  in  order  to  be  a  valid  representation  of  the  total  field. 


The  distorted-wave  perturbative  solutions  are  derived  with  the  assumption 
that  information  regarding  the  gross  (i.e.  purely  depth-dependent)  sound  speed 
structure  of  the  ocean  waveguide  is  known.  Using  this  knowledge,  a  background 
field,  P0(k0,r,rniz,zn),  which  satisfies 


|v2  +  klnl{z)\^P0{k0,r,ro,z,z0)  =  - 


2nr 


(2) 


can  be  calculated.  Since  n0(z)  is  known,  the  P0  field  can  be  thought  of  as  an 
incident  field  in  this  inhomogeneous  environment.  Because  the  oceanic  sound  speed 
structure  has  an  inherent  depth  variability,  but  depends  weakly  on  range,  the 
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background,  n0(z),  is  chosen  to  be  strictly  depth-dependent,  but  is  ’close*  (in  a 
sense  specified  below)  to  the  true  index  of  refraction,  n(r,z),  for  all  ranges  and 
depths.  By  perturbing  around  the  incident  field  solution,  approximate  solutions  to 
the  total  field  can  be  formed. 

2.3  The  Distorted- Wave  Born  (DWB)  Solution 

The  DWB  solution  is  formed  by  expanding  the  total  field,  F(k0,r,r0,z,z0),  and 
the  index  of  refraction,  n2(r,z),  in  powers  of  a  small  parameter,  e.  The  total  field, 
which  satisfies  Eq.  (1),  is  written  as 

P(Jc0,r,r0,z,z0)  =  P0  +  eP1  -h  •  •  •  (3) 

and  the  square  of  the  index  of  refraction  as 

n2(r,z)  =  n2(z)  +  enx(r,z)  +  •  •  •  (4) 

Substituting  Eqs.  (3)  and  (4)  into  Eq.  (1),  equating  powers  of  e  and  retaining 
terms  only  up  to  0(e)  results  in  two  equations.  The  0(e° )  equation  is  Eq.  (2). 
and  the  0(e )  equation  is 

V2  +  A:2n2(z)|p i(k0,r,r0,z,z0)  =  -klnl(r,z)P0  .  (5) 

Because  Eq.  (5)  is  the  Helmholtz  equation  with  a  source  term,  and  our  Green’s 
function  is  the  incident  field,  P„,  the  solution  to  Eq.  (5)  is  formed  by  integrating 
over  all  source  points.  The  field  solution  is  independent  of  theta  due  to  cylindrical 
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symmetry.  Performing  the  theta  integration  we  obtain 


P  i  {kot r,  r0yz,z0) 


oo  n 

=  2 vk20J  J  P0{k0tr,r\z,z')nl(r\z,)P0(k0tr',r0,z' ,z0)r‘dr'dz' 


(6) 


0  0 

for  the  first  term  in  the  perturbation  expansion.  Therefore,  to  first  order  in  e,  the 
total  field  is  given  by 


P(k0,r,rotz,zo)  ~  P[B)  {ko,ryr0,z,zo)  =  P0  +  eP\  ,  (7) 

where  P[B)  {k0,r,ro,z,zo)  is  defined  as  the  first  DWB  approximation  to  the  field 
(e.g.  Butkov  1968).  Because  Pj  is  commonly  called  the  scattered  field,  Eq.  (7)  is 
rewritten  as 


P{B)[ko,r,rp,zfz0)  =  P0  +  Pt . 


(8) 


where  c  has  now  been  set  equal  to  1. 

The  validity  of  the  approximation  is  based  on  the  assumption  that  the  total 
field  is  close  to  the  incident  field.  Equivalently,  |  P,/P  0  I  must  be  much  less  than 
1;  i.e.  the  scattered  field  is  small  compared  to  the  incident  field.  The  accuracy  of 
the  approximation  depends  on  the  size  of  the  perturbation,  the  propagation 
distance  in  the  perturbed  medium  and  the  source  frequency.  We  discuss  this 


further  below. 
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2.4  The  Distorted-Wave  Rytov  (DWR)  Solution 


The  first  DWR  approximation  is  simply  related  to  the  first  DWB  and  can 
therefore  be  derived  from  Eq.  (8)  (Fiddy  1986).  Dividing  both  sides  of  Eq.  (8)  by 
P0  and  taking  the  natural  logarithm  gives 


In 


Pl*> 


=  In 


1  + 


(9) 


or  since  ln(l+o:)~a  when  a«l,  we  can  define  the  first  DWR  approximation,  P[R ^ 
by 


In 


|>l*) 

Po 


(10) 


Exponentiating  both  sides  and  multiplying  by  P0  gives  the  first  DWR 
approximation 

P[RHK^,r0,z,z0)  ~  Po**'  •  (H) 

where  4>i  =  PJP  0  is  the  first  term  in  a  formal  phase  expansion, 

$  =  <f>\  +  <i>2  +  •  •  •  (12) 

and  4>i  is  O(t),  <£2  is  0(e2),  etc.  The  full  DWR  field  is  derived  in  Appendix  C. 

Because  we  are  mainly  interested  in  the  behavior  of  the  first  order 
approximations,  the  terms  DWB  and  DWR  will  denote  the  first  order 
approximations  unless  otherwise  indicated. 
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2.5  Discussion  of  Errors 


As  is  evident  from  Eqs.  (7)  and  (11)  and  the  derivation  of  the  DWR  from  the 
DWB,  the  first  order  solutions  axe  closely  related.  However,  as  derived  by  Keller 
(1969)  for  the  constant  background  solutions,  and  by  Belykin  (1985)  for  the 
distorted-wave  solutions,  the  first  order  relative  error  in  the  two  approximations  is 
quite  different.  By  relating  terms  in  the  series  in  Eqs.  (3)  and  (12),  Beylkin 
establishes  the  relative  DWB  error 


P  -  P[° ) 


=  0 


<f>2  ~  ~<£l2 


(13) 


while  the  relative  DWR  error  is 


P  -  P[*> 


=  o  j^J 


(14) 


Although  both  expressions  are  order  e2,  the  DWB  error  has  an  extra  term, 
<f>x2.  Since  (as  we  shall  see  later)  both  <f>i  and  <f>2  depend  on  range  to  the  first 
power,  the  DWB  error  term  is  0(e2r2),  while  the  DWR  is  0(c2r).  Because  of  this, 
the  amplitude  of  the  DWB  field  will  grow  like  r  as  we  propagate  in  range.  The 
DWR  is  expected  to  give  a  better  approximation  to  the  transmitted  field. 

We  must  keep  in  mind,  however,  that  these  are  only  first  order  errors  and,  in 
addition,  they  are  based  on  the  assumption  that  the  scattered  field  is  much  less 
than  the  incident  field.  As  previously  discussed,  the  DWB  requires  that 
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|  Pt/P  o  |  «1  in  order  for  the  approximation  to  be  valid.  But,  although  the 
DWR  was  derived  from  the  DWB  using  this  assumption,  the  latter  is  not  a 
necessary  condition  for  the  validity  of  the  DWR. 

However,  the  DWR  approximation  has  additional  terms  which  are  not  evident 
in  the  simple  derivation  of  the  DWR  from  the  DWB.  In  Appendix  C,  where  the 
full  DWR  approximation  is  derived,  the  assumption  is  made  that  the  integral  of 
[V$]2,  where  $  is  the  full  DWR  phase,  is  of  lower  order.  This  term  will  be  of  0(e2) 
only  if  the  phase  gradient  is  small;  we  can  have  no  sharp  discontinuities  in  the 
complex  phase  term.  This  makes  the  behavior  of  the  DWR  difficult  to  predict 
because  the  above  error  estimate  is  locally  inaccurate  if  the  ratio  P,/ P  0  is 
changing  rapidly. 

2.6  Explicit  Solutions 

The  solution  to  Eq.  (2)  can  be  written  in  several  ways.  In  this  thesis,  the 
incident  field  is  written  as  a  sum  of  discrete  modes.  The  solution  is  accurate  if  we 
are  in  the  far  field;  i.e.  the  continuous  spectrum  is  negligible  and  only  the  discrete 
modes  contribute  significantly  to  the  field.  The  modal  solution  for  an  arbitrary 
source  point,  derived  in  Appendix  A,  is 

Poi.^-oiriroyz\zo)  ~  ^  ^  (^wr)  J0{Cn  ro )  r>ro  (l5.a) 

n 
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Po{k(nryroi2,Zo)  —  ^  ^  i  £n)  0(^1  £n)  ^  (^nro),^o(^nr)  r<^ro  (15. b) 

n 

where  the  sum  over  n  includes  all  the  discrete  modes.  The  radial  wavenumber 
eigenvalue  is  £n  and  Bi2^((nr)  is  the  outgoing  wave  Hankel  function  of  zeroth  order 
with  harmonic  time  convention,  tlwl.  The  DWB  solution  is  written  as 

p[B)  {k0,r,r0,z,z0) 

00  k 

=  p0  +  2xk2J  j  P0(k0,r,r' ,z,z')ni(r\z')P0(k0,r' ,r0,z' >z0)r’ dr' dz' . 

0  O 

Substituting  Eqs.  (15)  into  (16),  and  setting  r0  =  0,  we  obtain 

,2  00  h 

PlB)(*„,r,0,z.z.)  =  P.  -  —■  f  f 

r  «n 

*Xl^%^M^,em)M2)(emr')rii(r',z')r'dr'da' 

m 

l2  r  h 

-  g  /  /  £#  ^,e„)^,e„)^o2)(Cnr)^(enr') 

0  o  n 

•J]V'(^.U^o,U^o2)(emr')ni(r',r')r'dr'dr'  .  (17) 


The  range  integration  can  be  done  analytically  if  the  true  profile  is  strictly 
depth-dependent  and  we  assume  this  here;  i.e.  we  set  nx{r' ,z')  =  nx[z').  The 
range  integrals,  which  we  define  as  I(n,m),  are  evaluated  for  two  cases;  Cm  =*= 
and  Cm  —  Cn  (see  Appendix  D). 
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If  in  Eq.  (17),  then  we  write  the  solution  symbolically  as 


P{B)  (fco,r,0 ,z,z0)  =  P0  +  Ptn  . 


(18. a) 


Pin ,  the  scattered  field  term  for  is  given  by 


ik^  r  ^ 

P»n{K>r,r0,z,z0)  =  -  f  ^rp{z' ,(n)ip(z,U)ni[z') 

o  n 


. «l!)(£„r)  ,  . 

nz  >Znmzo,tn) - - - dz  ■ 

Cn 


(18-b) 


where  H[ ^  (£„r)  is  the  first  order  Hankel  function.  If  Cm  ^  £n  then 


P[B){k0,r,Otz,z0)  ~  P0  +  P, 


snm 


(18.c) 


and  the  scattered  field  defined  as  Ptnm  is 


tic*  rh 

Psnm(^o>rtro>Z,Z0)  ~~7  J  ^  jtf>(z  £n)  ^1  (z  ) 


o  n 


V"\  ,/>+  \  1/  \  f^oHtmr)  ,  , 

'Zjnz  ,tmmzo,tm) - = - -dz 


m 


e  *  _  c  * 

tm  <»n 


(18.d) 


The  DWR,  formed  from  the  DWB  as  previously  described,  is 


P[R]  {K,r,0,z,zo)  =  P0exp 


sn 


Po 


(I9.a) 


if  tm  —  £n-  Similarly,  for  =*=  £„,  we  obtain 


P[R)  {ko,r,0,z,zo)  =  P 0exp 


snm 


(19-b) 
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In  the  next  chapter,  the  above  equations  will  be  evaluated  for  several  constant 
velocity  profiles  (i.e.  the  depth  dependent  solutions  arc  found  explicitly)  and  the 

behavior  of  the  Born  and  Rytov  in  these  simple  cases  is  investigated.  Although  not 
representative  of  the  velocity  structure  in  an  ocean  waveguide,  analyzing  these 
cases  gives  an  indication  of  how  the  DWB  and  DWR  approximations  will  behave 
in  a  complex  environment. 
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3.  CONSTANT  BACKGROUND  TEST  CASES 

3.1  Introduction 

In  this  chapter,  the  behavior  of  the  Born  and  Rytov  approximations  in  two 
constant  velocity  environments  is  discussed.  Although  a  constant  velocity 
waveguide  is  not  representative  of  a  deep  ocean,  the  solutions  and  error  estimates 
can  be  found  analytically.  Studying  the  behavior  of  the  approximations  using 
simple  velocity  models  often  gives  an  intuitive  feel  for  how  the  solutions  will 
behave  in  a  variable  velocity  waveguide. 

In  the  first  example,  the  true  and  guess  profile  are  chosen  to  be  constant  over 
the  entire  waveguide.  The  total  field  is  calculated  using  the  true  profile  and  the 
incident  field  using  the  guess  profile.  The  Born  and  Rytov  solutions  to  the  total 
field  are  calculated  and  compared  with  the  true  solution.  This  case  indicates  how 
the  approximate  solutions  will  behave  if  the  incident  or  guess  profile  is  different 
from  our  true  profile  over  the  width  of  the  waveguide,  or  equivalently,  dictates  how 
closely  our  guess  profile  must  match  the  true  when  the  profile  perturbation  extends 
over  the  full  waveguide. 

The  second  model  investigated  is  a  two-layer  waveguide.  The  true  profile  has 
two  constant  velocity  layers  and  the  first  guess  to  this  profile  is  a  constant  velocity 
waveguide.  Although,  again  not  representative  of  an  ocean  waveguide,  this  case 


provides  information  on  how  the  Bom  and  Rytov  approximations  will  behave  when 
the  difference  in  the  true  and  guess  fields  arises  from  a  small  velocity  perturbation 
(confined  to  a  portion  of  the  waveguide)  about  a  known  background.  The  goal  in 
this  thesis  is  to  recover  small  perturbations  about  a  depth-dependent  background 
profile,  and  approximate  estimates  of  the  magnitude  and  spatial  extent  of  allowable 
perturbations  for  the  variable  velocity  case  can  be  obtained  from  this  example. 

In  both  of  these  cases,  the  solutions  can  be  found  analytically  and  so  simple 
error  estimates  can  be  written  explicitly.  In  addition,  the  constant  velocity  cases 
serve  as  a  check  on  the  accuracy  of  the  numerical  calculations. 

3.2  Constant  Velocity  Waveguide 

3.2.1  Derivation  of  the  Solutions 

The  total  field  in  a  constant  velocity  waveguide,  n(z)  =  ne,  is  given  by 

P(ko,r,0,z,zo)  =  jsm\nz0smXnzH<^{^nr) ,  (20) 

n 

where  the  denote  the  horizontal  wavenumber  eigenvalues  for  the  total  field  and 
are  defined  by  =  ifc^n2  —  A„2.  The  reference  wavenumber  is  k0  =  —  and  the 

co 

vertical  wavenumbers  are  the  A„.  The  sine  eigenfunctions  arise  from  a  particular 
choice  of  boundary  conditions.  The  pressure  field  is  chosen  to  be  zero  at  the 
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surface  and  the  z-derivative  (since  we  have  a  flat  bottom)  to  be  zero  at  the  ocean 
bottom.  The  factor  (2 /h)  1 /2  is  the  normalization  constant  for  this  set  of 
orthonormal  eigenfunctions. 

The  initial  guess  field,  for  the  constant  velocity  waveguide  with  index  of 
refraction,  n0  is  given  by 


Poi.k0,  Zq)  ~  ^  sinAnz0  sinAnz  H\)  )  {^nf) , 


where  the  horizontal  wavenumber  eigenvalues  for  the  background  solution  axe 


given  by  the  £n  which  are  defined  by  £n2  =  Jk2n$  —  An2.  As  is  evident  from  Eqs. 


(20)  and  (21),  the  eigenfunctions  are  identical.  In  a  constant  velocity  waveguide, 
the  only  difference  in  the  two  solutions  are  the  values  of  the  horizontal 
wavenumbers. 

From  Eq.  (17),  the  Born  solution  in  a  constant  velocity  medium  is  given  by 


P[D]  {K,r,Q,z,z0)  = 


Ik2  f  h 

JE  sinA„zsinAnz'  H[2)  (£„r)  J0(£„r ') 

^  o  o  n 


*2  sinA™2'  sin*™2o  M2)  {tmr')r‘dr'dz' 


k20nx-n 


oo  h 


/  /  E  sinA„*sinA„z'  H[2]  {Znr')J0{Znr) 


r  o  n 


’2  sin^m«'sinAmz0  H$(Zmr')r’dr'dz'., 


where  nlt  the  constant  velocity  perturbation,  is  the  difference,  n2  -  n2,  in  the 
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squares  of  the  true  and  background  refractive  indices. 

Because  the  perturbation  is  constant,  both  the  range  and  depth  integrals  can 
be  evaluated  explicitly.  Evaluating  the  range  integrals  (see  Appendix  D)  and 
following  the  notation  of  Eqs.  (18),  the  Bom  scattered  field  for  is 


•  I  2  » 

ik0riyr  r  ™ 

Psn  (k0,r,ro,z,z0)  =  — — ^  J  £jsinXnZsm 


2  A nz'  sinA, 


M2)  (e,r) 

e„ 


dz'  . 


(23.a) 


For  =*=  £m,  we  obtain 


P»nm{koiri  roizizo) 
l2  h 


ik0n  1 


f  V~\  •  ,  r-.  •  ,  •  #o2^(£mr) 

J  Sln^nzsm^n2'  2_j  sinAmz'  sinAmz0 - - - —  dz' 


o  n 


m 


C 


(23.b) 


Because  the  eigenfunctions  are  orthogonal  over  the  width  of  the  waveguide  and 
there  is  no  mode  coupling  due  to  the  absence  of  a  depth-dependent  perturbation, 
Eq.  (23. b)  is  always  zero.  Evaluating  the  depth  integral  in  Eq.  (23.a),  results  in 
the  constant  velocity  scattered  field  solution 

.  .  .  , 

sinA„zsinAnz0 - - - 

tn 

Therefore,  the  Bom  and  Rytov  solutions  axe  formed  by  substituting  Eq.  (24)  into 
Eqs.  (18.a)  and  (19. a),  respectively. 

The  first  order  errors  in  the  two  solutions  are  discussed  in  the  next  section. 


Ptn{koirirotzizo) 


Ah 


Because  the  errors  have  the  form  of  a  ratio  of  sums,  explicit  expressions  are  found 
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for  a  single  mode  only.  However,  as  will  be  seen  in  the  numerical  examples,  the 
single  mode  error  estimates  closely  approximate  the  error  for  the  multi-modal 
solutions. 


3.2.2  Error  Analysis 


As  discussed  in  section  2.5,  the  first  order  errors  can  be  derived  in  terms  of 
the  first  two  Rytov  phase  terms,  ^  and  <^>2.  The  first  phase  term  in  a  constant 
velocity  waveguide  is  just 


ik?0ni  r 

4>l{ko,r,0,z,zo)  =  -  ~  ~  2_j  sinAnzsinAng0 


AhP, 


H\2][jnr) 

in 


(25) 


The  second  phase  term,  fa,  is  derived  from  the  expression  for  the  full  Rytov 
solution  (see  Eq.  C.7).  Assuming  that  the  perturbation  difference,  nl5  is  known 
exactly,  then  to  first  order  fa  is  given  by 


oo  h 


fa(ko,r,0,z,zo)  =  jrf  f  P0(K^r\z,z') 


0  0 


\^<h{koyr\0,z'  ,z0)  P0{ko,r\r0,z',z0)r'dr'dz' 


(26) 


therefore,  the  gradient  of  4>i  and  the  incident  field,  P0(k0,r,Ofz,z0),  are  needed  to 
calculate  <^2. 
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Because  the  phase  terms  have  the  form  of  a  ratio  of  sums,  a  simple  error 
estimate  can  not  be  established  for  all  modes.  However,  if  we  choose  to  look  at  the 
errors  for  a  single  mode  only,  Po(kotr,0,z,zo)  can  be  written  explicitly  as 

Poi(ko,r,0,z,zo)  =  —smXiZsmXiZo^i^r)  (27) 

where  P0i[k0,r,Otz,z0)  denotes  the  field  expression  for  the  first  mode,  Aj  is  the 
corresponding  vertical  wavenumber  while  6  is  the  corresponding  horizontal 
wavenumber.  This  notation  will  be  used  throughout  in  deriving  the  subsequent 
single  mode  expressions. 


Using  Eq.  (27),  the  first  phase  term  can  now  be  written 


$11  (^01  t),z,z0)  — 


f. 


(28) 


Using  the  asymptotic  forms  of  the  Hankel  functions 


M2)(e  ir)  = 


*tl r 
2 

7r^r 


^  -.-(6  r-f) 


(29) 


and  canceling  terms,  <frn(ko,r,0,z,zo)  is  given  by 


ikjnvr 

26 


$li  {ko,r,0,z,zo)  = 


(30) 


so  the  gradient  of  <f>u  is  just 
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il?0n  i 

'V4>n{k0,r,Q,z,z0)  =  r  . 

-‘Cl 


(31) 


Substituting  Eq.  (31)  into  Eq.  (26)  and  evaluating  the  range  integral  for  a  single 
mode,  <^>2i  becomes 


<t>-n{K,r,Q,z,z0)  =  - 


k\nxr  M2)(fiO 


*0*1 

26 


2  A 


J  sin2Ai z‘ dz' .  (32) 


Performing  the  depth  integration  and  cancelling  like  factors,  the  second  phase  term 
becomes 


irk\n\ 

86: 


<M*<Mr>°,‘*>2o)  =  - 


(33) 


and  we  can  now  write  the  first  order  error  estimates. 

From  Eq.  (13),  the  first  order  Born  error  expression  for  a  single  mode  is 


p-p[B)  0 

r2k\n\ 

+  0 

irk*  n\ 

p 

862 

8£i3 

(34) 


and  from  Eq.  (14),  the  first  order  Rytov  error  estimate  is 


P-P[R ) 


=  0 


irk\n\ 

86: 


(35) 


As  discussed  in  section  2.5,  the  Born  error  has  a  term  that  is  of  the  order  r2nf, 
while  the  Rytov  error  is  order  rnf.  The  growth  in  the  Born  approximation  to  the 
transmitted  field  is  well-known  (e.g.  Keller,  1969)  and  the  Rytov  solution  is 
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expected  to  give  a  better  approximation  to  the  true  field  solution. 

Although  the  first  order  error  estimates  suggest  in  general  how  the  two 
solutions  will  behave,  a  modal  analysis  of  the  fields  in  the  next  section,  yields  the 
errors  in  the  first  order  solutions  exactly  for  a  single  mode.  The  errors  are  seen 
through  the  projection  of  the  mode  coefficient. 


3.2.3  Modal  Decomposition  of  the  Bom  and  Rytov  Fields 


Because  the  background  solution  is  a  modal  representation  of  the  field,  the 
Born  and  Rytov  representations  to  the  total  field  can  also  be  considered  to  be 
composed  of  a  sum  of  normal  modes.  In  general,  the  true  mode  functions  would  be 
the  unknowns,  but  in  a  constant  velocity  waveguide  the  functional  form  of  the 
eigenfunctions  for  the  guess  profile  are  identical  to  those  of  the  true  profile.  This 
permits  a  straightforward  evaluation  of  the  errors  in  the  Born  and  Rytov  solutions. 

From  Eq.  (E.5),  the  Born  coefficient  for  a  single  mode  is 


and  Eq.  (E.6)  for  the  Rytov  coefficient  is 


(  ' 

2 

1/2 

7T  Za 

irk20n , 

2  fi 

ClR  = 

h 

\  4 

Sm(~2A_)eXP 

(37) 


Clearly,  if 


n ! 
2fi 


«1,  then  the  two  coefficients  are  the  same  to  first  order.  If, 
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however,  this  quantity  is  not  small,  the  Bom  exhibits  an  algebraic  amplitude 
growth,  while  the  Rytov  error  is  oscillatory.  Therefore,  although  the  Born  error 
estimates  in  section  3.2.2  are  valid  for  any  profile  perturbation,  the  Rytov  error 
estimate  is  correct  only  if  the  scattered  field  is  small  with  respect  to  the  incident 
field. 

In  the  next  section,  the  numerical  results  for  a  constant  velocity  waveguide 
are  presented  and  the  analytical  results  of  this  section  will  be  verified. 

3.3  Constant  Velocity  Waveguide  Example 

The  parameters  for  the  constant  velocity  waveguide  example  are  given  in 
Table  1.  The  true  and  background  profiles  are  illustrated  in  Figure  7.  The 
transmission  loss  curves  for  the  profiles  are  illustrated  in  Figure  8.  The 

I  P  I 

transmission  loss  is  defined  as  —  20log(-j-1 — L— ),  where  P  is  the  complex  field  and 

I  PrtJ  I 

Prej  is  the  field  at  1  meter  from  the  source. 

Because  the  perturbation  is  small  and  the  eigenfunctions  are  identical  in  the 
constant  velocity  case,  the  true  and  background  curves  nearly  overlay  each  other, 
except  perhaps  at  the  ''nulls”  where  the  modes  have  interfered  destructively 
creating  a  minimum  amplitude.  The  amplitude  behavior  of  the  Born  and  Rytov 
solutions  can  be  qualitatively  predicted  by  studying  Figure  9;  the  ratios  j  PJP  0  | 
vs.  range.  The  ratios  match  until  approximately  5  km;  this  is  the  regime  in  which 
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Table  1.  Input  Parameters:  Test  Case  1 


INPUT  PARAMETERS:  TEST  CASE  1 

Range 

1.-50.  km 

Source  Depth 

500.  m 

Receiver  Deptii 

500.  m 

Source  Frequency 

25.  Hz 

Depth  (m) 

Sound  Speed  (m/s) 

True  profile  0. 

1500. 

1500. 

1500. 

Guess  profile  0. 

1501. 

1500. 

1501. 

Maximum  Perturbation 

1.  m/s 

Perturbation  Width 

1500.  m 

the  ratio  is  less  than  1/4.  Beyond  this  range  the  Born  ratio  continues  to  grow 
while  the  Rytov  peaks  at  a  ratio  of  two  and  then  begins  to  decrease.  Although  the 
behavior  of  the  Born  solution  is  easily  predicted,  determining  the  behavior  of  the 
Rytov  is  more  difficult. 

The  Rytov  ratio  does  not  exhibit  algebraic  growth  except  at  small  ranges 
where  both  solutions  are  within  their  range  of  validity.  From  the  modal 
decomposition  in  section  3.2.3,  we  could  assume  little  or  no  amplitude  error. 
However,  neglecting  the  gradient  of  the  ratio  of  the  scattered  to  incident  field  in 


T-3790 


40 


Figure  7. 


SOUND  SPEED  vs.  DEPTH 


Constant  velocity  profiles.  The  true  sound  speed  is  1500  m/s.  The 
guess  sound  speed  in  1501  m/s. 


TRANSMISSION  LOSS  vs.  RANGE 


True  and  background  transmission  loss  curves  for  test  case  1.  The 
curves  axe  generated  using  the  profiles  in  Figure  7.  Because  the 
eigenfunctions  have  the  same  functional  form,  the  curves  are 
indistinguishable. 


Figure  8. 
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Ps  /PQ  vs.  RANGE 

O 


RANGE  (km) 

Figure  9.  Ratio  of  the  scattered  field  to  the  incident  field  for  the  Bom  and 
Rytov  solutions  for  test  case  1.  The  solutions  are  generated  using 
the  background  profile  in  Figure  7. 


this  decomposition  will,  in  some  cases,  cause  local  errors  in  the  Rytov  field.  In  this 
example,  we  might  be  inclined  to  predict  a  local  error  in  the  field  solution  at  the 
spike  near  SO  km.  However,  estimating  the  value  of  the  local  gradient  near  this 
point,  results  in  a  value  less  than  1.  Therefore,  the  Rytov  solution  is  predicted  to 
everywhere  match  the  true  solution. 

The  Bom  solution  to  the  field  and  the  true  solution  are  illustrated  in  Figure 
10.  The  growth  of  the  Bom  solution  with  range  is  as  predicted  from  the  error 
analysis  and  the  modal  decomposition.  The  Rytov  solution  to  the  field  and  the 
true  solution  are  illustrated  in  Figure  11.  As  expected  from  the  results  in  the 
previous  section  and  the  above  discussion,  the  Rytov  exhibits  no  amplitude  growth, 
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Figure  10. 


Figure  11. 


TRANSMISSION  LOSS  vs.  RANGE 


O 


True  and  Born  transmission  loss  curves  for  test  case  1.  The  growth 
of  the  Born  solution  with  range  is  clearly  seen. 


TRANSMISSION  LOSS  vs.  RANGE 


True  and  Rytov  transmission  loss  curves  for  test  case  1.  Unlike  the 
Born  solution,  the  Rytov  is  unaffected  by  the  extent  of  the 
perturbation. 
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and  matches  the  true  profile  almost  exactly. 


The  behavior  of  the  approximations  is  clearly  demonstrated  by  viewing  the 
relative  error  in  the  transmission  loss  of  the  background,  the  Born  and  Rytov 
solutions  (see  Figure  12)  with  respect  to  the  true  transmission  loss.  The  Rytov 
solution  is  the  best  fit;  the  small  variations  of  the  guess  from  the  true  have  been 


smoothed.  However,  for  large  r,  the  relative  error  in  the  Born  transmission  loss 


2log(nl£lr) 

grows  approximately  like  - ; — -  as  we  propagate  in  range.  This  estimate  is 

logfir 


obtained  from  the  single  mode  solution.  For  a  single  mode,  the  amplitude  behavior 

with  range  for  the  true  transmission  loss  is  roughly  log[l/(£ ,  r)1/2  ]  and  for  the 


RELATIVE  ERROR  vs.  RANGE 


1 .0  25.5  50  0 


RANGE  (km) 

Figure  12.  Relative  error  in  the  background,  Born  and  Rytov  transmission  loss 
curves  for  test  case  1.  The  Born  error  is  approximately  17%  at  50 
km.  The  error  in  the  Rytov  solution  is  approximately  zero  over  the 
entire  range  of  propagation. 


Born  transmission  loss  is  log[l/(£  if)1/2  (l+r2fc04n12/^  i2)1^2  ].  By  assuming 
k0~£i  and  canceling  terms,  we  obtain  the  above  expression  for  the  relative  error. 

To  further  confirm  the  results  of  the  modal  decomposition,  we  present  some 
numerical  examples.  Figures  13  and  14  illustrate  the  results  of  the  Born  and  Rytov 
modal  decompositions,  respectively,  for  a  single  mode.  The  amplitude  of  the  Born 
coefficient  and  the  real  part  of  the  Rytov  coefficient  confirm  the  validity  of  the 
expressions  for  the  mode  coefficients  in  the  previous  section.  The  amplitude  of  the 
Born  coefficients  grows  with  range  while  the  error  in  the  real  part  of  the  Rytov 
coefficients  oscillate  between  0  and  200  percent.  On  the  other  hand,  the  error  in 
the  real  part  of  the  Born  coefficients  remains  stable  about  zero  as  does  the 
amplitude  of  the  Rytov  coefficients. 

To  confirm  that  the  multi-mode  behavior  matches  that  of  the  single  mode,  a 
modal  decomposition  of  the  solutions  using  the  first  seven  modes  is  illustrated  in 
Figures  15  and  16  over  a  range  of  1-20  km.  So,  even  for  multiple  modes,  the  errors 
predicted  by  the  single  mode  decomposition  are  qualitatively  correct.  The 
accuracy  of  the  multi-mode  decompositions  is  verified  by  forming  the  transmission 
loss  using  the  recalculated  Born  and  Rytov  coefficients. 

The  real  part  of  the  Born  coefficients  axe  used  to  recalculate  the  Born 
transmission  loss.  The  resultant  curve,  labeled  'MDB'  is  illustrated  in  Figure  17 
along  with  the  true  transmission  loss  and  the  original  Born  solution.  As  expected, 
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o  RELATIVE  ERROR:  BORN  COEFFICIENT 


RANGE  (km) 


Figure  13.  Relative  error  in  the  amplitude  and  real  part  of  the  Born  mode 
coefficient  for  a  single  mode  in  test  case  1.  The  real  part  of  the  Born 
coefficient  is  identically  equal  to  the  true  coefficient. 


o  RELATIVE  ERROR:  RYTOV  COEFFICIENT 

o 


RANGE  (km) 


F(gure  14.  Relative  error  in  the  amplitude  and  real  part  of  the  Rytov  mode 
coefficient  for  a  single  mode  in  test  case  1.  The  amplitude  of  the 
Rytov  coefficient  is  identically  equal  to  the  true  coefficient. 
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Figure  15. 


Figure  16. 


o  RELATIVE  ERROR:  BORN  COEFFICIENTS 


Relative  error  in  the  amplitude  and  real  part  of  the  Born  mode 
coefficient  for  the  First  seven  modes  in  test  case  1. 


9  RELATIVE  ERROR:  RYTOV  COEFFICIENTS 


Relative  error  in  the  amplitude  and  real  part  of  the  Rytov  mode 
coefficient  for  the  first  seven  modes  in  test  case  1. 


TRANSMISSION  LOSS  vs.  RANGE 


Figure  17.  True,  Born  and  the  reconstructed  Born  transmission  loss  curves 
generated  using  the  first  seven  modes  in  test  case  1.  The 
reconstructed  Born  solution  is  calculated  using  the  real  part  of  the 
projected  mode  coefficients. 


TRANSMISSION  LOSS  vs.  RANGE 


Figure  18.  True,  Rytov  and  the  reconstructed  Rytov  transmission  loss  curves 
generated  using  the  Fust  seven  modes  in  test  case  1.  The 
reconstructed  Rytov  solution  is  calculated  using  the  amplitude  of  the 
projected  mode  coefficients. 
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the  growth  of  the  Born  solution  with  range  has  disappeared.  To  recalculate  the 
Rytov  solution,  we  use  the  amplitude  of  the  mode  coefficients.  The  result, 
illustrated  in  Figure  18,  and  labeled  ”MDR'  exhibits  little  change  from  the  original 
Rytov  solution. 

In  reviewing  the  results  of  this  section,  we  conclude  that  although  the 
amplitude  of  the  Born  approximation  is  strongly  affected  by  a  shift  in  the  sound 
speed  over  the  full  waveguide,  the  Rytov  approximation  is  insensitive  to  a 
perturbation  of  this  magnitude  and  extent.  In  the  next  section,  we  study  several 
two-layer  waveguides  in  order  to  determine  the  magnitude  and  extent  of 
perturbations  recoverable  about  a  known  background  profile. 

3.4  Two-Layer  Waveguide 

3.4.1  Derivation  of  the  Solutions 

The  two-layer  model  used  in  this  section  is  illustrated  in  Figure  19.  The 
sound  speed  in  the  upper  layer  is  defined  to  be  Cj,  while  the  sound  speed  in  the 
lower  layer  is  defined  to  be  c0,  to  match  the  constant  velocity  case.  The  reference 
wavenumbers  are  kl  and  k0,  respectively.  The  depth  of  the  layer  is  Zj.  The  total 


field  for  this  geometry  is 
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P(k,r,0,z,zo) 


i  \  \  ,  cos\on{h-z0)cos\oll{h-Zi) 

=  - — “ — 3 - smAlnz//J,2'(fnr)  0  <z<zx 


4  "  sinAlnzj 

n 

F(A^r,0,z,zo) 

=  A*COsX«n(h-Z0)cos\on{h-z) 


zi<z<h  (38) 


where  An  is  given  by 
K  -  2AlflAonsin2AtBr, 

•  |AMco8aAM(A-2i)  {z,Alft  ~  sinAlflz,cosAuziJ 

+  AInsin2Ajnz,  [(A-Zi)A,„  +  sinA^A-ZiJcosA^A-zj)]  J  (39) 

Again,  the  fB  are  the  horizontal  eigenvalues  for  the  total  Field,  and  are  defined  by 


TWO —LAYER  WAVEGUIDE  GEOMETRY 


z  =  0 


*1 


z  =  h 


Cl 


Figure  19.  Two-layer  geometry  used  in  test  cases  2  and  3. 
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fn2  =  klnj  -  A?„  0 <z<z1\  =  k20n20  -  A2n  zx<z<h. 

The  initial  guess  field,  calculated  using  a  constant  velocity  waveguide,  with  index 
of  refraction,  n0 ,  is  given  in  Eq.  (21).  Note,  that  although  the  eigenvalues  in  the 
lower  layer  in  the  total  field  are  given  in  terms  of  k0,  they  are  not  identical  with 
those  in  the  constant  velocity  case.  The  eigenvalues  for  the  background  profile  are 
still  given  by  £fl2  =  /c2n2  -  A2. 

The  Born  solution  is  given  by  Eqs.  (18)  where  the  integration  is  now  over  the 
support  of  the  perturbation  in  depth,  from  0  to  zx.  Because  the  background 
eigenfunctions  are  not  orthogonal  over  the  perturbation  extent  in  depth,  Eq.  (I8.d) 
is  non-zero  in  this  example.  The  depth  integral  can  still  be  evaluated  in  this  case. 
The  scattered  Born  solution  for  is  given  by 


F>sn[k:oiriroi!'izo)  ~ 


ik20  n  i  r 
2/i2 


Y  ~  u7sin2A”2‘ 


•  sinA„zsinA„z0 


In 


(40. a) 


For  the  solution  is 


f3tnm{koirir<nZ,Z<l) 


ik20nxr 


J]sinAnzsinA„z0 


^(Ur) 

Zm'-Zn2 


XmcosXnz1  sinA^Zj  -  AnCosA^ZjsinA^Z] 
A  2  -  A  2 

An  Am 


(40.b) 


Analogous  to  the  constant  velocity  example,  the  Born  solution  is  formed  by 


substituting  the  scattered  field  solutions,  Psn  and  Pgnm,  into  Eqs.  (18.a)  and  (18.c), 
respectively;  the  Rytov  solution  is  formed  by  substituting  Pgn  and  Psnm  into  Eqs. 
(I9.a)  and  (l9.b),  respectively. 

The  errors  in  the  two-layer  waveguide  solutions  are  discussed  in  the  next 
section.  As  in  the  previous  example,  the  explicit  error  expressions  will  be  derived 
for  a  single  mode  only  and  these  simple  estimates  will  be  compared  to  the  results 
from  the  multi-modal  solutions. 


3.4.2  Error  Analysis 


Following  the  analysis  in  section  (3.2.2),  we  will  derive  the  first  order  errors 
for  the  case  of  a  single  mode.  From  Eq.  (40.a),  the  first  phase  term  in  a  two-layer 
waveguide  is 


r, 0,2, - - - J  sm(— )dz.  (41) 


Using  P0  from  Eq.  (21),  the  asymptotic  form  of  the  Hankel  functions  (Eq.  (29)), 
and  evaluating  the  integral,  we  obtain 


<t>n{K,r>0,z,zo)  = 


iklrrij 

2£i 


.  **  l  , 

4,  S1"(  — ' 


7T 


(42) 


The  gradient  of  4>n  is  just 
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Wii(*o»r,0,*,zo)  = 


*1  n  t 
2£i 


A 


•  ( irZl  \ 
sin(~T) 

7T 


so,  from  Eq.  (26),  fci  simply  becomes 


jkfrn  i 


<^2l(^o>r>0,i,Zo)  =  - 


•  f*Zi\ 

«.  sm(T> 


From  Eqs.  (42)  and  (44),  the  first  order  relative  Bom  error  is  given  by 
P-P[B] 


=  o 


ik04rnj  /c04  r2  nf  ^ 


8&3  8^2 

and  the  first  order  relative  Rytov  error  is 


nzi  ^  M 


«,  sin(T') 


7T 


P-P[R)  „ 

- 1 

. 

o** 

I* 

-t 

MM 

k  ^1 

2‘ 

p  -° 

8«13 

A  7T 

<  J 

(43) 


(44) 


(45) 


(46) 


The  behavior  of  the  first  order  relative  errors  in  range  thus  depends  critically  on 
the  width  of  the  layer  in  depth.  If  zj  is  small  compared  to  A,  the  sine  can  be 
expanded  around  Zi/h  =  0  and  the  errors  now  look  like 

?H*.  )61 

(47) 


p-p[b)  Q 

<KArn\  k04r2  n\ 

z 1 1 

6' 

P 

sei3  1  h2  . 

1 h  \ 

for  the  Born  error  and 


T-3790 


53 


p-p is>  c 

ik0 4  mj 
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6 

p 

Hi 3 

h  1 

l  J 

(48) 


for  the  Rytov  error.  Therefore,  if  the  depth  of  the  layer  is  small  with  respect  to 
the  depth  of  the  waveguide,  the  two-layer  waveguide  errors  will  be  decreased  from 
the  constant  velocity  errors  by  a  factor  ( z^/h )  6.  From  these  results,  it  would  be 
expected  that  the  errors  in  both  solutions  would  decrease  for  any  layer  width  less 
than  the  full  extent  of  the  waveguide. 


3.4.3  Modal  Decomposition  of  the  Bom  and  Rytov  Fields 

In  contrast  with  the  constant  velocity  example,  the  true  and  background 
mode  functions  are  not  identical  in  this  case.  The  background  eigenfunctions  axe 
sines  but  the  true  mode  functions  are  more  complicated  due  to  the  presence  of  the 
layer.  Correspondingly,  the  expressions  for  the  mode  coefficients  are  more  difficult 
to  derive  and  interpret. 

From  Eq.  (E.8),  the  Born  coefficient  for  the  single  mode  in  the  two-layer 
waveguide  is 


and  from  Eq.  (E.10),  the  Rytov  is 
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nzn 


Cir  =  2Asin(-^-)exp 


ik\rnx 

2hix 


iL 

h 


•  (*ZI\ 
sm(~) 

7T 


(50) 


where  A  is  given  in  Eq.  (E.9).  As  in  the  constant  velocity  case,  the  Rytov  mode 
coefficient  reduces  to  the  Born  if 


iklrni 

2htx 

h  7T 

«  1 


and  we  truncate  the  series  to  first  order. 

Superficially,  the  modal  decompositions  seem  to  suggest  the  same  type  of 
behavior  as  the  first  order  errors.  However,  as  will  be  seen  in  the  next  section,  the 
quantity  A  plays  a  significant  role  in  determining  the  validity  of  the  solutions. 


3.5  Two-Layer  Waveguide  Examples 

In  this  section,  we  discuss  two  velocity  models;  the  first  has  a  100  m  layer  and 
the  second  a  400  m  layer.  The  parameters  for  the  first  model  are  given  in  Table  2. 

The  true  and  background  profiles  are  illustrated  in  Figure  20  and  the 
transmission  loss  curves  in  Figure  21.  Again,  the  guess  overlays  the  true  profile 
almost  exactly  due  to  the  smallness  of  the  perturbation  and  also  its  minimal 
extent. 

We  will  first  qualitatively  predict  the  behavior  of  the  Born  and  Rytov  solutions 
by  studying  the  ratio  )  Pt/P  0  J  vs.  range  in  Figure  22.  Both  of  the 
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Table  2.  Input  Parameters:  Test  Case  2 


INPUT  PARAMETERS: 

TEST  CASE  2 

Range 

1.-50.  km 

Source  Depth 

500.  m 

Receiver  Depth 

500.  m 

Source  Frequency 

25.  Hz 

Depth  (m) 

Sound  Speed  (m/s) 

True  profile  0. 

1500. 

100. 

1500. 

100. 

1501. 

1500. 

1501. 

Guess  profile  0. 

1501. 

1500. 

1501. 

Maximum  Perturbation 

1.  m/s 

Perturbation  Width 

100.  m 

ratios  are  less  than  1/4  (except  at  the  spikes)  and  they  match  in  amplitude  except 
at  the  spikes  which  occur  at  the  nulls  in  the  amplitude  of  the  incident  field.  This  is 
a  case  where  we  axe  in  the  range  of  validity  of  both  solutions;  the  ratio  of  the 
scattered  to  the  incident  field  is  small.  Therefore,  the  solutions  are  expected  to 
generate  similar  results. 

The  Born  and  the  true  solution  are  illustrated  in  Figure  23  and  the  Rytov  and 
the  true  solution  are  illustrated  in  Figure  24.  As  expected,  the  approximate 
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SOUND  SPEED  vs.  DEPTH 


1498  1500  1502 


Figure  20.  The  true  and  background  sound  speed  profiles  used  in  test  case  2. 

The  background  profile  is  constant  over  the  entire  waveguide  at  1501 
m/s.  The  true  profile  has  a  100  meter  layer  in  which  the  sound 
speed  is  1500  m/s.  Below  the  layer,  the  true  and  background  souno 
speeds  are  identical. 


TRANSMISSION  LOSS  vs.  RANGE 


O 


Figure  21. 


The  true  and  background  transmission  loss  curves  for  test  case  2. 
Because  the  perturbation  is  small  (l  m/s)  and  extends  only  100 
meters  in  depth,  the  two  curves  are  nearly  identical. 
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Figure  22. 


Figure  23. 


Ps  /PQ  vs.  RANGE 


tatio  of  the  scattered  field  to  the  incident  field  for  the  Born  and 
Rytov  solutions  in  the  two-layer  waveguide  of  test  case  2.  Because 
the  trend  of  the  ratio  is  less  than  1/4  over  the  range  of  propagation, 
both  solutions  are  expected  to  yield  accurate  results. 


TRANSMISSION  LOSS  vs.  RANGE 


The  true  and  Born  transmission  loss  curves  for  test  case  2.  The  Born 
solution  matches  the  true  except  at  the  nulls  in  the  field. 
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TRANSMISSION  LOSS  vs.  RANGE 


i 
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RANGE  (km) 

Figure  24.  The  true  and  Rytov  transmission  loss  curves  for  test  case  2.  The 
Rytov  solution  matches  the  true  except  at  the  nulls  in  the  field. 


RELATIVE  ERROR  vs.  RANGE 


1.0  25.5  50.0 


RANGE  (km) 


Figure  25.  Relative  error  in  the  background,  Born  and  Rytov  transmission  loss 

curves  for  test  case  2.  The  error  in  the  the  solutions  is  always  less 
than  10%. 
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solutions  overlay  the  true  solution,  except  at  the  nulls,  as  is  illustrated  in  Figure 
25,  by  the  relative  errors. 

At  this  point,  we  could  further  confirm  the  above  results  by  proving  the  errors 
to  be  minimal  in  the  modal  decomposition,  but  instead  we  will  just  conclude  that 
for  a  1  m/s  perturbation  extending  over  7%  of  the  waveguide  the  Born  and  Rytov 
are  valid  approximations  over  a  range  of  1  to  50  km. 

The  next  case  discussed  is  that  of  a  layer  extending  over  28%  (400  m)  of  the 
waveguide.  The  input  parameters  for  this  example  are  given  in  Table  3. 

The  true  and  background  profiles  are  illustrated  in  Figure  26  and  the 
transmission  loss  curves  are  illustrated  in  Figure  27.  There  is  a  visible  difference 
in  the  true  and  guess  curves.  From  test  case  1,  we  know  that  when  the 
perturbation  is  1  m/s  over  the  full  waveguide,  the  transmission  loss  curves  are 
basically  identical.  The  reason  for  the  variation  in  test  case  3  is  due  to  the  change 
in  the  eigenfunctions.  As  the  width  of  the  layer  increases  to  a  certain  depth,  the 
difference  in  the  true  and  background  mode  functions  increases  correspondingly. 
As  the  layer  width  extends  to  the  full  waveguide,  the  mode  shapes  return  to  those 
found  in  a  constant  velocity  environment. 

Figure  28  demonstrates  the  change  in  the  first  mode  when  the  layer  width  is 
increased  in  200  m  increments.  The  first  curve  is  the  mode  shape  in  a  constant 
velocity  waveguide;  the  second  curve  is  the  mode  shape  when  a  200  m  layer  is 


T-3790 


60 


Table  3.  Input  Parameters:  Test  Case  3 


INPUT  PARAMETERS:  TEST  CASE  3 

Range 

1.-50.  km 

Source  Depth 

500.  m 

Receiver  Depth 

500.  m 

Source  Frequency 

25.  Hz 

Depth  (m) 

Sound  Speed  (m/s) 

True  profile  0. 

1500. 

400. 

1500. 

400. 

1501. 

1500. 

1501. 

Guess  profile  0. 

1501. 

1500. 

1501. 

Maximum  Perturbation 

1.  m/s 

Perturbation  Width 

400.  m 

present,  the  third  curve  is  the  mode  shape  when  a  400  m  layer  is  present  and  so 
on. 

To  estimate  the  percentage  change  in  the  mode  functions,  we  look  at  the 
change  in  the  mode  coefficient  over  the  waveguide  as  the  layer  width  is 
systematically  increased.  From  Figure  29,  we  can  estimate  that  a  400  meter  layer 
with  a  1  m/s  change  in  the  sound  speed,  results  in  approximately  a  14%  change  in 
the  mode  coefficient  as  compared  to  the  constant  velocity  case.  It  is  evident,  then, 
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SOUND  SPEED  vs  DEPTH 
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f  igure  26.  The  true  and  background  sound  speed  profiles  used  in  test  case  3. 

The  background  profile  is  constant  over  the  entire  waveguide  at  1501 
m/s.  The  true  profile  has  a  400  meter  layer  in  which  the  sound 
speed  is  1500  m/s.  Below  the  layer,  the  true  and  background  sound 
speeds  are  identical. 


TRANSMISSION  LOSS  vs.  RANGE 


o 


I 


1.0  25.5  50  0 

RANGE  (km) 

Figure  27.  The  true  and  background  transmission  loss  curves  for  test  case  3. 

Although  the  perturbation  is  small  (1  m/s),  it  extends  over  28%  of 
the  waveguide;  therefore  the  guess  transmission  loss  differs  from  the 

true. 
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MODE  SHAPE  vs.  LAYER  WIDTH 


o 


Figure  28.  Illustrated  is  the  variation  of  the  lowest  order  mode  as  a  function  of 
layer  depth.  The  curve  on  the  far  left  represents  the  mode  shape  in 
a  constant  velocity  waveguide  and  the  second  curve  from  the  left  is 
the  mode  shape  when  a  200  m  layer  is  present,  etc.  The  curve  on 
the  far  right  is  the  shape  of  the  lowest  order  mode  when  a  1.4  km 
layer  is  present.  When  the  layer  extends  to  1.5  km,  the  shape 
returns  to  that  in  a  constant  velocity  waveguide. 

MODE  COEFFICIENT  vs.  LAYER  WIDTH 


Figure  29.  The  change  in  the  mode  coefficient  for  the  lewest-order  mode  as  a 
function  of  the  layer  width.  In  a  constant  velocity  waveguide,  the 
numerical  value  is  —0.018.  The  value  of  the  coefficient  decreases  to 
~.005  at  1  km,  and  at  1.5  km  returns  to  the  original  value  of  0.018. 
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that  a  small  change  in  the  sound  speed  over  28%  of  the  waveguide  causes  a  large 
change  in  the  value  of  the  eigenfunctions  and  a  corresponding  change  in  the  true 
and  background  transmission  loss  curves. 

To  predict  the  behavior  of  the  Born  and  Rytov  solutions,  we  again  study  the 
ratio  of  the  scattered  to  the  incident  field  for  both  approximations.  From  Figure 
30  and  a  closeup  in  Figure  31,  we  r  that  the  ratio  in  both  cases  has  exceeded  1/4 
by  a  range  of  approximately  13  km.  Therefore,  the  Born  solution  is  expected  to 
grow  for  ranges  greater  than  13  km.  The  behavior  of  the  Rytov  solution  is  difficult 
to  predict.  Although  the  Rytov  ratio  also  exceeds  1/4,  we  saw  in  the  constant 
velocity  case  that  a  large  ratio  did  not  affect  the  validity  of  the  Rytov  solution. 

The  Born  and  true  transmission  loss  plots  are  illustrated  in  Figure  32  and  the 
Rytov  and  true  in  Figure  33.  As  expected,  there  is  a  slight  growth  in  the  Born 
solution.  The  Rytov  solution  is  somewhat  more  accurate  (see  the  relative  error 
versus  range  in  Figure  34),  but  again  does  not  match  the  true  solution;  in  addition, 
a  spike  in  the  solution  appears  at  about  5  km  in  the  same  location  as  one  of  the 
spikes  in  the  |  Pt/P  0  |  of  Figure  30. 

We  now  see  the  effect  that  a  14%  change  in  the  true  and  background 
eigenfunctions  has  on  the  accuracy  of  the  Born  and  Rytov  solutions.  We  also  see 
the  failure  of  the  error  predictions  to  account  for  this  type  of  change.  Only  the 

factor  \/~A  in  the  modal  decompositions  gives  an  indication  that  the  solutions  will 
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Figure  32. 


Figure  33. 


TRANSMISSION  LOSS  vs.  RANGE 


The  true  and  Bom  transmission  loss  curves  for  test  case  3.  Although 
the  Born  reconstructs  the  general  shape  of  the  curve,  the  growth  of 
the  solution  is  apparent  for  ranges  greater  than  ~13  km. 
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The  true  and  Rytov  transmission  loss  curves  for  test  case  3. 
Although  the  Rytov  approximates  the  true  solution  better  than  the 
Born,  amplitude  errors  are  apparent  in  both  the  nulls  and  peaks  in 
the  field. 
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RELATIVE  ERROR  vs.  range 
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Figure  34.  Relative  error  in  the  background,  Born  and  Rytov  transmission  loss 
curves  for  test  case  3.  Although  the  Rytov  has  the  least  error, 
neither  of  the  perturbative  solutions  were  able  to  duplicate  the  true 
transmission  loss  curve. 


become  invalid  for  layer  widths  less  than  the  full  waveguide. 

Although  difficult  to  analyze  quantitatively,  the  significance  of  \[a  can  be 
explained  in  terms  of  the  constant  velocity  versus  the  two-layer  solution.  From 
Eq.  (21),  we  know  the  true  and  background  solutions  in  a  constant  velocity 
waveguide  are  constructed  from  an  identical  set  of  eigenfunctions;  only  the 
eigenvalues  differ.  In  the  two-layer  waveguide,  the  eigenfunctions  of  the  true 
solution  vary  (see  Eq.  (38))  depending  on  the  receiver  location  in  the  waveguide 
with  respect  to  the  layer.  The  former  case  has  a  simple  normalization  factor  of 
v(2//i) ,  while  from  Eq.  (39)  the  latter  has  the  expression,  .  From  the  modal 
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decomposition  of  the  Bom  and  Rytov  approximations  (Appendix  E) ,  the  term  V  A 
occurs  which  contains,  in  addition  to  the  factor  V~A,  other  terms  resulting  from  the 
projection  of  the  coefficients  over  a  depth  0  to  zj  where  the  constant  velocity 
eigenfunctions  are  not  orthogonal.  Therefore,  for  the  Born  and  Rytov  solutions  to 

be  valid  representations  of  the  true  field,  the  term  \Ta  must  give  the  right  mode 

shapes  for  the  two-layer  waveguide.  This  is  impossible  because  \f=A  is  constant  for 
a  given  layer  depth  and  there  still  remains  only  a  single  set  of  eigenfunctions 
predicting  the  true  behavior  over  the  entire  waveguide.  As  a  result,  the  Born  and 
Rytov  solutions  are  unable  to  approximate  the  true  solution  if  the  true  and  guess 
eigenfunctions  differ  by  more  than  about  7%. 

The  results  of  this  section  will  now  be  used  to  predict  the  behavior,  in  the 


next  chapter,  of  the  perturbation  expansions  in  a  depth-dependent  environment. 


T-3790 


68 


4.  DEPTH-DEPENDENT  BACKGROUND  TEST  CASES 

4.1  Introduction 

Although  a  constant  background  solution  provides  useful  error  estimates,  it  is 
not  viable  zeroth  order  solution  in  the  deep  ocean.  For  this  reason,  we  now  turn  to 
the  distorted-wave  perturbative  approximations  in  which  the  background  profile 
depends  on  depth.  Several  depth-dependent  velocity  profiles  are  studied  and  the 
results  compared  to  those  predicted  from  the  constant  velocity  analysis.  In  all 
cases,  the  trend  of  the  sound  speed  profile  is  assumed  to  be  known  exactly;  we  look 
at  the  behavior  of  the  DWB  and  DWR  solutions  when  a  perturbation  is  placed 
upon  the  background  solution.  Unlike  the  previous  chapter,  the  results  in  this 
chapter  cannot  be  derived  explicitly  even  for  a  single  mode.  The  reason  is  that  the 
eigenfunctions  are  known  only  numerically.  The  DWB  and  DWR  representations 
of  the  total  field  are  obtained  from  Eqs.  (7)  and  (10),  respectively. 

The  background  profile  in  the  first  two  examples  is  the  weakly  refractive 
NORDA  2A  test  case  (NORDA  Parabolic  Workshop  1984).  Two  different  size 
perturbations  (referred  to  as  test  cases  4  and  5)  are  placed  on  this  profile  and  we 
study  the  ability  of  the  DW  approximations  to  construct  the  true  solution. 

In  the  second  set  of  examples,  the  profile  is  strongly  refractive.  Two  different 
test  cases  (referred  to  as  test  cases  6  and  7)  are  presented.  The  perturbation  is  the 
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same  is  both  cases  but  the  source  position  and,  therefore,  the  illumination  in  the 
waveguide  is  varied. 


4.2  NORDA  2 A  Examples 


The  parameters  used  in  test  case  4  are  given  in  Table  4. 


Table  4.  Input  Parameters:  Test  Case  4 


INPUT  PARAMETERS: 

TEST  CASE  4 

Range 

1.-50.  km 

Source  Depth 

500.  m 

Receiver  Depth 

500.  m 

Source  Frequency 

25.  Hz 

Depth  (m) 

Sound  Speed  (m/s) 

True  profile  0. 

1500.0 

300. 

1506.0 

400. 

1507.5 

500. 

1510.0 

1000. 

1520.0 

1500. 

1563.0 

Guess  profile  0. 

1500.0 

1000. 

1520.0 

1500. 

1563.0 

Maximum  Perturbation 

0.5  m/s 

Perturbation  Width 

200.0  m 
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The  NORDA  2A  bilinear  profile,  illustrated  in  Figure  35,  is  the  background 
profile  in  this  example.  The  true  profile  is  constructed  by  placing  a  small 
perturbation  of  limited  extent  on  top  of  the  guess.  The  profile  perturbation  is 
illustrated  in  Figure  36;  over  the  rest  of  the  profile,  the  true  and  background  sound 
speed  values  are  identical. 

The  transmission  loss  curves  for  the  perturbed  and  unperturbed  profiles  are 
illustrated  in  Figure  37.  Even  though  the  profile  perturbation  is  only  .5  m/s  and 
only  extends  over  13%  of  the  waveguide,  there  are  visible  differences  between  the 
two  curves.  However,  from  the  results  of  section  3.5,  we  expect  the  DWB  and 
DWR  to  give  a  good  approximation  to  the  true  solution. 


SOUND  SPEED  vs.  DEPTH 


Figure  35.  NORDA  2A  bilinear  background  profile  used  in  test  cases  4  and  5. 
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Figure  36. 


Figure  37. 


SOUND  SPEED  vs.  DEPTH 


1506  1508  1510 


The  profile  perturbation  used  in  test  case  4.  The  perturbation 
extends  from  300  to  500  m  in  depth  and  the  maximum  variation 
from  the  background  is  .5  m/s.  In  the  remainder  of  the  waveguide, 
the  true  sound  speed  matches  the  background. 
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TRANSMISSION  LOSS  vs.  RANGE 
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True  and  background  transmission  loss  curves  for  test  case  4. 
Although  the  maximum  difference  in  sound  speed  is  only  0.5  m/s, 
the  perturbation  extent  is  ~13%  of  the  waveguide;  therefore  the 
curves  differ  as  we  propagate  in  range. 
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The  ratios  of  the  scattered  to  the  incident  field  amplitudes  as  a  function  of 
range  are  illustrated  in  Figure  38.  In  order  to  clearly  illustrate  the  trend  of  the 
ratios,  the  same  curves  are  shown  in  Figure  39  with  an  expanded  vertical  scale. 
From  Figure  38,  we  would  predict  the  possibility  of  a  large  local  error  in  the  DWR 
solution  due  to  the  spike  near  40  km.  Figure  39,  on  the  other  hand,  gives  an 
indication  of  the  global  behavior  of  the  solutions,  particularly  the  DWB.  Because 
the  trend  of  the  ratio  does  not  exceed  1/4,  we  would  expect  to  obtain  an  accurate 
solution  from  the  DWB  as  well  as  the  DWR.  In  fact,  for  this  test  case,  the 
expansions  are  approximately  equal  for  ranges  less  than  50  km. 

Figures  40  and  41  illustrate  the  transmission  loss  results  from  the  DWB  and 
DWR  solutions.  The  DWB  exhibits  minimal  growth  for  ranges  greater  than  about 
40  km,  while  the  DWR  solution  contains  a  spike  at  the  null  at  the  same  range. 
But  as  illustrated  by  the  relative  errors  in  Figure  42,  both  solutions  give  a  good 
approximation  to  the  true  transmission  loss. 
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Figure  38. 


Figure  39. 


Ps  /PQ  vs.  RANGE 


Ratio  of  the  scattered  field  to  the  incident  field  for  the  DWB  and 
DWR  solutions  for  test  case  4.  The  spike  at  ~40  km  suggests  the 
possibility  of  a  local  error  in  the  solutions. 


o  Ps  /P0  vs.  RANGE 
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Expanded  version  of  Figure  38.  The  trend  of  the  ratios  is  less  than 
1/4  and  the  DWB  and  DWR  are  predicted  to  give  accurate  results. 
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Figure  40. 


Figure  41. 


TRANSMISSION  LOSS  vs.  RANGE 


RANGE  (km) 

True  and  DWB  transmission  loss  curves  for  test  case  4.  As  predicted, 
the  DWB  is  a  good  approximation  to  the  true  solution. 
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True  and  DWR  transmission  loss  curves  for  test  case  4.  As 
predicted,  the  DWR  is  a  good  approximation  to  the  true  solution 
except  at  the  null  near  40  km. 
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RELATIVE  ERROR  vs.  RANGE 


1.0  25  5  50.0 


RANGE  (km) 


Figure  42.  Relative  error  in  the  background.  Bom  and  Rytov  transmission  loss 
curves  for  test  case  4. 


In  the  next  example,  the  size  of  the  perturbation  is  increased  to  1  m/s  but  the 
extent  of  the  perturbation  remains  the  same.  The  input  parameters  for  test  case  5 
are  given  in  Table  5.  The  background  solution  is  again  the  NORDA  2A  profile. 

The  profile  perturbation  is  illustrated  in  Figure  43.  The  transmission  loss 
curves  for  the  true  and  background  profiles  are  illustrated  in  Figure  44.  As 
expected,  the  difference  in  the  true  and  guess  curves  is  greater  than  in  the  previous 
example.  The  guess  and  the  true  solutions  differ  both  in  magnitude  and  in  the 
position  of  the  nulls. 

The  scattered  to  incident  field  ratios  are  illustrated  in  Figure  45  and  a  closeup 
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Table  5.  Input  Parameters:  Test  Case  5 


INPUT  PARAMETERS:  TEST  CASE  5 

Range 

1.-50.  km 

Source  Depth 

500.  m 

Receiver  Depth 

500.  m 

Source  Frequency 

25.  Hz 

Depth  (m) 

Sound  Speed  (m/s) 

True  profile 

0. 

1500. 

300. 

1506. 

400. 

1507. 

500. 

1510. 

1000. 

1520. 

1500. 

1563. 

Guess  profile 

0. 

1500. 

1000. 

1520. 

1500. 

1563. 

Maximum 

Perturbation 

1.  m/s 

Perturbation  Width 

200.  m 

of  the  same  curves  is  illustrated  in  Figure  46.  In  this  case,  the  magnitude  of  the 
spike  (although  not  indicated  by  the  plot  scale)  in  Figure  45  near  40  km  is  over 
twice  that  in  the  previous  example  and  we  would  expect  a  larger  local  error.  From 
Figure  46,  we  would  expect  the  DWB  solution  to  fail  at  ranges  greater  than  25  km 
when  the  ratio  exceeds  1/4.  Again,  the  transmission  loss  curves  illustrated  in 
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SOUND  SPEED  vs.  DEPTH 

1506  1508  1510 


Figure  43.  Profile  perturbation  for  test  case  5.  The  maximum  variation  from 
the  background  is  now  1  m/s  and  the  extent  is  ~13%  of  the 
waveguide. 
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TRANSMISSION  LOSS  vs.  RANGE 


Figure  44. 


True  and  background  transmission  loss  curves  for  test  case  5.  Note 
the  difference  in  the  two  curves  as  compared  to  Figure  37. 
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Ps  /P0  vs.  RANGE 


Figure  45.  Ratio  of  the  scattered  field  to  the  incident  Field  for  the  DWB  and 
DWR  solutions  for  test  case  5.  As  in  Figure  38,  the  spike  at  ~40  km 
suggests  the  possibility  of  a  local  error  the  solutions. 


Ps  /PQ  vs.  RANGE 


Figure  46.  Expanded  version  of  Figure  45.  The  DWB  solution  is  expected  to  be 
in  error  past  25  km  when  the  ratio  exceeds  1/4. 
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Figures  47  and  48  verify  the  qualitative  predictions.  Although,  the  DWB  solution 
has  shifted  from  the  guess  to  indicate  correctly  the  position  of  the  nulls,  the  growth 
in  the  solution  becomes  evident  past  25  km.  The  DWR,  on  the  other  hand, 
contains  a  spike  at  the  null  near  40  km  and,  in  addition,  incorrectly  positions  the 
nulls.  The  amplitude  behavior  is,  however,  closer  to  the  true  solution  than  the 
DWB.  The  relative  error  in  the  transmission  loss  of  the  guess,  DWB  and  DWR 
solutions  is  illustrated  in  Figure  49. 

This  is  exactly  the  behavior  predicted  from  the  two-layer  case  in  section  3.5. 
As  the  change  in  the  true  and  background  mode  functions  increases,  the  ability  of 
the  DWB  and  DWR  to  construct  the  true  solution  decreases. 


TRANSMISSION  LOSS  vs.  RANGE 


RANGE  (km) 

Figure  47.  True  and  DWB  transmission  loss  curves  for  test  case  5.  Although 
the  shape  of  the  curve  is  correct,  the  growth  in  the  DWB  solution  for 
ranges  greater  than  25  km  is  evident. 
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TRANSMISSION  LOSS  vs.  RANGE 
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Figure  48.  True  and  DWR  transmission  loss  curves  for  test  case  5.  The  error  in 
the  DWR  solution  near  ~40  km  has  increased  (compare  Figure  41) 
and  there  are  errors  in  the  amplitude.  In  addition,  the  DWR  fails  to 
correctly  position  the  nulls  of  the  solution. 

RELATIVE  ERROR  vs.  RANGE 


Figure  49.  Relative  error  in  the  background,  DWB  and  DWR  transmission  loss 
curves  for  test  case  5.  The  vertical  scale  has  been  expanded  to 
differentiate  the  errors. 
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We  now  turn  to  some  deep  ocean  examples  in  which  the  sound  energy 
propagates  mainly  through  refraction. 

4.3  SOFAR  Channel  Examples 

As  discussed  in  Chapter  1,  the  velocity  structure  of  the  SOFAR  channel  is 
characterized  by  a  sound  speed  minimum  (known  as  the  channel  axis)  at  about  1 
km  with  the  sound  speed  increasing  both  above  and  below  the  axis.  In  this  section, 
we  study  the  behavior  of  the  DWB  and  DWR  solutions  when  the  perturbations  are 
placed  about  a  bilinear  background  sound  speed  profile.  In  particular,  we 
investigate  the  effect  of  varying  the  source  position  with  respect  to  the  channel 
axis. 

The  parameters  for  the  first  example  in  this  section,  test  case  6,  are  given  in 
Table  6.  The  source  is  placed  at  500  m,  while  the  channel  axis  is  at  1  km. 
Because  the  profile  is  strongly  refractive,  the  entire  waveguide  is  not  insonified. 
Instead  the  energy  in  the  waveguide  travels  in  a  relatively  narrow  beam. 

The  background  profile  is  illustrated  in  Figure  50  and  the  perturbation  about 
the  background  is  illustrated  in  Figure  51.  The  true  and  background  transmissions 
loss  curves,  illustrated  in  Figure  52,  give  an  indication  of  the  insonified  and 
shadow  zones  in  the  waveguide.  Near  the  source  the  energy  is  a  maximum,  but  as 
we  move  out  in  range  to  approximately  13  km  there  is  a  broad  low  energy  region 
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Table  6.  Input  Parameters:  Test  Case  6 


INPUT  PARAMETERS:  TEST  CASE  6 

R 

ange 

1.-50.  km 

Source  Depth 

500.  m 

Receiver  Depth 

500.  m 

Source 

Frequency 

25.  Hz 

Depth  (m) 

Sound  Speed  (m/s) 

True  profile 

0. 

1520. 

700. 

1478. 

800. 

1473. 

900. 

1466. 

1000. 

1460. 

4000. 

1525. 

Guess  profile 

0. 

1520. 

1000. 

1460. 

4000. 

1525. 

Maximum 

Perturbation 

1.  m/s 

Perturn 

ition  Width 

200.  m 

which  extends  to  approximately  30  km.  Past  30  km,  the  region  around  the 
receiver  depth  is  again  illuminated. 

The  ratios  of  he  scattered  to  incident  field,  illustrated  in  Figure  53  with  an 
closeup  in  Figure  '>1,  indicate  that  the  maximum  error  is  expected  to  occur  in  the 
shadow  zone;  in  p.  .  ticular,  in  the  transition  from  low  to  high  energy  near  a  range 
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SOUNO  SPEED  vs.  DEPTH 
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Figure  50.  Background  profile  used  in  test  cases  6  and  7. 


SOUND  SPEED  vs.  DEPTH 
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Figure  51.  Profile  perturbation  for  test  cases  6  and  7.  The  maximum  variation 
from  the  background  is  1  m/s  and  the  extent  is  ~5%  of  the 
waveguide. 
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TRANSMISSION  LOSS  vs.  RANGE 


Figure  52.  True  and  background  transmission  loss  curves  for  test  case  6.  Note 
the  broad  region  of  low  intensity. 


of  30  km.  However,  except  for  isolated  spikes  the  ratio  is  always  less  than  1/4  and, 
therefore  the  perturbative  solutions  are  expected  to  yield  accurate  results. 

The  DWB  and  DWR  transmission  loss  curves  axe  illustrated  in  Figures  55  and 
56.  The  vertical  scales  on  both  plots  have  been  expanded  to  include  the  spike 
occurring  near  29  km.  The  DWB  and  DWR  solutions  axe  nearly  identical  as 
predicted  by  the  small  |  PJP  „  |  ratios.  Both  solutions  match  the  true  solution  in 
the  regions  of  high  intensity,  but  fail  to  give  the  correct  values  in  the  shadow  zone 
and  in  fact  exhibit  a  greater  error  than  the  guess  near  the  deep  null.  A  plot  of  the 
relative  errors  (Figure  57)  in  the  two  solutions  substantiates  this  conclusion. 
Although  qualitative  predictions  suggest  that  both  solutions  axe  within  their  range 
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R.  /PQ  vs.  RANGE 


Figure  53.  Ratio  of  the  scattered  Field  to  the  incident  field  for  the  DWB  and 
DWR  solutions  for  test  case  6. 


/P0  vs.  RANGE 


Figure  54. 


Expanded  version  of  Figure  53.  Note  the  maximum  error  occurs  in 
the  region  of  low  intensity. 
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TRANSMISSION  LOSS  vs.  RANGE 


Figure  55.  True  and  DWB  transmission  loss  curves  for  test  case  6.  The  vertical 
scale  has  been  expanded  to  include  the  spike  near  29  km. 


TRANSMISSION  LOSS  vs.  RANGE 


Figure  56.  True  and  DWR  transmission  loss  curves  for  test  case  6.  The  vertical 
scale  has  been  expanded  to  include  the  spike  near  29  km. 
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RELATIVE  ERROR  vs.  RANGE 
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Figure  57.  Relative  error  in  the  background,  DWB  and  DWR  transmission  loss 
curves  for  test  case  6. 


of  validity,  this  test  case  indicates  that  the  DWB  and  DWR  are  invalid  over  broad 
regions  of  low  intensity  as  well  as  at  the  isolated  nulls  in  the  Held. 

In  the  last  example,  the  parameters  are  the  same  as  in  test  case  6  except  the 
source  and  receiver  are  placed  on  the  channel  axis  at  I  km.  The  true  and  guess 
transmission  loss  plots  are  illustrated  in  Figure  58.  The  shadow  zone  has 
disappeared  as  expected  and  the  background  solution  nearly  overlays  the  true. 

As  illustrated  in  Figure  59,  the  DWB  and  DWR  scattered  to  incident  field 
ratios  are  nearly  identical.  A  closeup  of  the  plots  in  Figure  60  demonstrates  the 
ratio  of  the  trends  to  be  less  than  1/4.  Although  the  spike  at  approximately  42  km 
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TRANSMISSION  LOSS  vs.  RANGE 
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Figure  58.  True  and  background  transmission  loss  curves  for  test  case  7. 

Because  the  source  is  placed  on  the  channel  axis,  the  shadow  zone 
has  disappeared. 


indicates  the  possibility  of  a  local  error,  the  DWB  and  DWR  solutions  are 
predicted  to  give  an  accurate  representation  of  the  true  transmission  loss.  This  is 
confirmed  by  the  transmission  loss  curves  in  Figures  61  and  62,  and  the  relative 
errors  in  Figure  63. 

The  results  of  this  chapter  verify  that  the  behavior  of  the  DWB  and  DWR 
solutions  can  be  qualitatively  predicted  from  the  quantitative  single  mode  results  of 
chapter  3.  If  the  eigenfunctions  do  not  change  by  more  than  approximately  7% 
and  the  receiver  is  well-illuminated,  the  perturbative  solutions  will  give  an  accurate 
representation  of  the  true  field  In  the  ocean  waveguide  to  a  range  of  50  km. 
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Figure  59.  Ratio  of  the  scattered  field  to  the  incident  field  for  the  DWB  and 
DWR  solutions  for  test  case  7. 


Ps  / PQ  vs.  RANGE 


Figure  60. 


Expanded  version  of  Figure  59.  The  ratio  is  always  less  than  1/4  and 
the  DWB  and  DWR  are  expected  to  give  accurate  results. 


T-3790 


90 


TRANSMISSION  LOSS  vs.  RANGE 
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Figure  61.  True  and  DWB  transmission  loss  curves  for  test  case  7.  Because  the 
waveguide  is  well-illuminated  and  the  perturbation  is  small,  the 
DWB  gives  a  good  approximation  to  the  true  solution. 


TRANSMISSION  LOSS  vs.  RANGE 
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Figure  62.  True  and  DWR  transmission  loss  curves  for  test  case  7.  Because  the 
waveguide  is  well-illuminated  and  the  perturbation  is  small,  the 
DWR  gives  a  good  approximation  to  the  true  solution. 
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Figure  63. 


RELATIVE  ERROR  vs.  RANGE 


Relative  error  in  the  background,  DWB  and  DWR  transmission  loss 
curves  for  test  case  7. 
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5.  INVERSION  RESULTS 

5.1  Introduction 

Inverting  for  a  depth-dependent  perturbation  in  the  ocean  waveguide  differs 
from  the  seismic  case  in  which  the  location  of  and  the  impedance  contrast  across 
sharp  boundaries  in  a  layered  earth  is  the  goal  (Cohen  and  Bleistein  1977;  Raz 
1981;  Beydoun  and  Tarantola  1988).  It  also  differs  from  the  problem  of  object 
scattering  for  which  the  definition  of  the  boundaries  of  a  compact  object  placed  in 
a  constant  background  are  desired  (Mueller,  Kaveh  and  Wade  1979;  Slaney,  Kak 
and  Larsen  1984;  Zapalowski,  Leeman  and  Fiddy  1985). 

In  the  ocean,  the  perturbations  about  the  inhomogeneous  background  are 
smoothly  varying  functions  of  depth  and  extend  over  the  entire  range  of 
propagation.  The  field  structure  in  the  deep  ocean  waveguide  is  characterized  by 
multiply  refracted  energy  as  opposed  to  reflected  energy.  Therefore,  if  we  think  of 
the  propagation  paths  in  terms  of  rays,  the  raypaths  can  cross  through  a 
perturbation  many  times  before  reaching  the  receiver.  Although  in  the  seismic 
case,  multiples  (energy  that  has  been  reflected  more  than  once  from  an  acoustic 
boundary)  are  inherent  in  the  data,  they  are  removed  by  processing  in  order  to 
display  the  primary  reflection;  energy  that  has  traveled  from  the  source  to  the 
layer  back  to  the  receiver  directly.  So  while  the  multiply  reflected  energy  is 
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removed  in  seismic  processing,  the  multiply  refracted  energy  in  the  ocean  is  the 
primary  arrival. 

The  applicability  of  Fourier  inversion  techniques  for  recovering  depth- 
dependent  perturbations  in  the  ocean  waveguide  is  investigated  in  this  chapter. 
Our  goal  is  to  relate  the  perturbation  as  a  function  of  depth  and  the  scattered  field 
as  a  function  of  frequency  as  a  Fourier  transform  pair  within  a  linear 
transformation.  As  will  be  shown,  the  DWR  and  DWB  scattered  field  data  are 
constructed  differently,  but  in  both  cases  the  transformation  kernel  is  obtained 
from  the  DWB  representation.  Therefore,  we  concentrate  on  the  attributes  of  the 
DWB  kernel;  the  behavior  of  the  amplitude  and  phase  derivative  with  respect  to 
depth  as  a  function  of  depth,  range  and  frequency.  Certain  criteria  must  be  met  in 
order  for  the  transform  relation  to  be  valid;  the  amplitude  must  be  a  slowly 
varying  function  of  frequency  and  the  phase  must  be  monotonic  in  depth  and  its 
frequency  dependence  known. 

In  the  following  sections,  we  discuss  the  inversion  formalism,  the  attributes  of 
the  DWB  kernel  and  the  effect  that  band-limiting  and  phase  averaging  have  on  the 
reconstruction  of  a  perturbation. 

5.2  The  Inversion  Algorithm 


The  DWB  scattered  field  representation  is  given  in  Eqs.  (18.b)  and  (18. d). 
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Because  the  solution  is  given  as  a  sum  over  all  modes,  explicit  expressions  for  the 
amplitude  and  phase  are  not  available.  If,  however,  we  perform  the  summation,  we 
can  write  the  scattered  field,  Pt,  as 


f 

Psi^o i 0, z, Zg'j  —  J  B(k0,r,z,z0,z  )t 


9(k„r,z,z0,z') 


n,(z')  dz‘  , 


(53) 


where  zx  and  z2  define  the  limits  of  the  perturbation  in  depth.  The  amplitude  of 
the  modal  sum  is  defined  as  B,  and  the  phase  as  9.  Be 16  is  defined  as  the 
transformation  kernel  relating  the  scattered  field  to  the  profile  perturbation. 

To  construct  the  transform  pair 


FT 

p,(K)  <=>  iiM*')), 


(54) 


we  generate  the  scattered  data  as  a  function  of  wavenumber  using  either  DWB  or 
DWR  approximations.  We  then  assume  that  the  amplitude  wavenumber 
dependence  is  known  and  can  be  absorbed  into  the  scattered  data  and  replace 
B(k0,r,z,z0,z’)  with  its  wavenumber  independent  equivalent,  B(r,z,z0,z');  the 
modified  scattered  field  is  written  as  Pt.  Next,  we  multiply  both  sides  of  Eq.  (53) 
by  the  phase  function,  e  integrate  over  wavenumber  to  obtain 


f  Pe(kotr~0,z,zo)e  'e^r^z-)dko 
—00 

z-l  00  _ 

=  J  nx  (z')B(r, z,z0,z‘)  J  e"(*- *-•*■*.•*')  "  dkg  dz‘  .  (55) 


-00 


Assuming  that  the  phase,  0,  depends  linearly  on  wavenumber,  9  =  k09  where  9  is 
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independent  of  k0  (or  if  not  a  transformation  has  been  made  to  account  for  the 
dependence),  and  the  phase  is  monotonic  in  z',  then  the  k0— space  integral  is 
identically  equal  to  a  delta  function.  The  argument  of  the  6-function  is  the 
difference  of  our  ^-independent  phase  terms,  O(z')  and  O(z').  Because  the  delta 
function  is  in  terms  of  the  phase  instead  of  depth,  we  make  a  transformation 
obtaining  the  integrand  as  a  function  of  the  phase.  As  a  result,  the  transform 
variables  are  k0  and  0(z').  Therefore,  we  write  Eq.  (55)  as 


w  .  _ 

f  Pt(k0)e-ikMz']dk0 

—OO 


r~ 


r 

Bz  | 


h,  nx{0{z))B(e{z))2*6\0{z')  -  6{z) 


I  dO/dz  '  |  z=- 


dO. 


(56) 


where  dOjdz  '  is  the  Jacobian  of  the  transformation.  For  simplicity,  we  have 
excluded  the  parameter  dependence  in  the  expression;  we  retain  only  the  transform 
variables  O(z')  and  k0  explicitly.  Evaluation  of  the  phase  integral  gives 


_  OO 


.,<*•))  = 

2nB(0(z')) 


P,{h)e 


dk„ 


(57) 


Thus,  we  have  obtained  a  linear  transform  relation  between  the  scattered  field  data 
as  a  function  of  wavenumber  and  the  perturbation  as  a  function  of  the  phase. 

To  form  the  DWB  scattered  field,  we  generate  the  true,  P,  and  background, 
P0,  field  representations  as  a  function  of  wavenumber  and  subtract  thereby 
obtaining 
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P[B](k0)  =  P(ko)-  Po(k0).  (58) 

The  DWR  scattered  field  is 

d*'(y  =  p-(*-)to7r^-  <59) 

The  DWB  and  DWR  representations  of  the  scattered  field  will  closely  approximate 
the  data  generated  from  Eq.  (53)  if  the  solutions  are  within  their  range  of  validity 
and  the  inversion  criteria  discussed  above  are  satisfied.  In  the  next  section,  we 
investigate  the  behavior  of  the  amplitude  and  the  phase  derivative  of  the  DWB 
kernel  for  two  sound  speed  profiles.  We  also  outline  the  processing  steps  necessary 
to  formulate  the  inversion  and  the  effects  that  these  processing  steps  have  on  the 
output  of  the  Fourier  inversion. 

5.3  Investigation  of  the  DWB  kernel 

The  first  case  investigated  is  the  NORDA  2A  profile  whose  parameters  where 
listed  previously  in  Table  1.  A  pixel  map  of  the  transmission  loss  over  a  range  of  1 
to  50  km  is  illustrated  in  Figure  64.  The  source  frequency  is  50  Hz.  Because  the 
sound  speed  minimum  is  located  at  the  ocean  surface,  there  are  multiple  surface 
reflections  interfering  with  the  refracted  energy  in  the  waveguide. 

As  indicated  in  section  5.1,  it  is  necessary  to  determine  the  behavior  of  the 
amplitude  and  phase  derivative  of  the  DWB  kernel  as  a  function  of  wavenumber. 
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Figure  65  illustrates  the  amplitude  of  the  DWB  kernel  while  Figure  66  illustrates 
the  phase  derivative  arising  from  the  transformation  from  z'  to  O(z'). 

Comparing  Figures  64  and  65,  we  see  that  they  contain  similar  but  shifted 
features.  The  shift  is  most  apparent  in  the  location  of  the  lower  turning  points  of 
the  field.  The  range  dependence  of  the  background  field  is 
=  J0(£r)  —  »T0(£r)  while  the  range  dependence  of  the  scattered  field  is 
iff)  (r)  =  Ji(fr)  -  iYiiCr).  In  the  far  field,  Ji(fr)~K0(^r)  and  J„(£r)  and 
K1(fr)  are  ir  out  of  phase  (Abramowitz  and  Stegun  1984). 

The  structure  in  Figure  66  is  essentially  the  instantaneous  vertical 
wavenumber,  kz,  analogous  to  the  instantaneous  frequency  in  complex  seismic  trace 
analysis  (Taner,  Koehler  and  Sheriff  1979).  The  yellows  and  reds  denote  positive 
values,  black  is  the  zero  transition  and  the  blues  are  negative  values.  The 
association  of  the  phase  derivative  with  kz  is  verified  by  the  zero  transition  regions 
in  Figures  65  and  66.  For  instance,  at  the  deep  turning  points  of  the  amplitude, 
the  phase  derivative  is  zero;  i.e.  the  propagation  is  horizontal.  Many  other  zero 
crossings  are  evident  in  the  figure  and  each  corresponds  to  kz  =  0. 

Recalling  that  the  phase  needs  to  be  monotonic  for  a  valid  inverse,  we  see  that 
at  least  globally  this  criterion  cannot  be  met.  However,  if  we  can  find  a  region 
where  the  phase  is  locally  monotonic,  and  the  rest  of  the  inversion  criteria  are  met, 
then  we  can  attempt  a  perturbation  reconstruction  over  a  limited  range  in  depth. 
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Figure  64. 


Transmission  loss  for  the  NORDA  2A  profile  illustrated  in  Figure  35. 
The  source  is  placed  at  a  depth  of  .5  km  and  the  source  frequency  is 
50  Hz. 


T-3790 


99 


Figure  65.  Amplitude  of  the  DWB  kernel  for  the  NORDA  2A  proli'.<  .  I  lie 
source  is  placed  at  a  depth  of  .5  km  and  the  source  frequency  is  50 
Hz.  Note  the  position  of  the  lower  turning  points  as  compared  to 
Figure  M. 
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Figure  66.  Derivative  of  the  phase  with  respect  to  depth  for  the  NORDA  2 A 
profile.  The  source  is  placed  at  a  depth  of  .5  km  and  the  source 
frequency  is  50  Hz.  The  complicated  pattern  of  zero  crossings  is  due 
to  the  interference  of  reflected  and  refracted  energy. 
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The  next  step  is  to  check  the  frequency  dependence  of  the  amplitude  and  phase 
derivative  or  equivalently,  the  phase,  at  a  fixed  range  (the  phase  is  constructed 
from  the  phase  derivative  to  avoid  phase  unwrapping).  The  amplitude  was  found 
to  have  a  \(k0  approximate  dependence  and  so  are  multiplied  by  k0.  The  phase 
derivative,  on  the  other  hand  goes  approximately  like  k0 ,  and  accordingly,  is 
divided  by  ka.  The  resultant  amplitude  and  phase  derivative  curves  are  illustrated 
in  Figures  67  and  68,  respectively,  for  10  frequency  values  ranging  from  5  to  50  Hz. 
Although  the  curves  in  each  plot  are  of  the  same  order  of  magnitude,  they  are  not 
identical.  Therefore,  only  frequency-averaged  versions  of  these  functions  are 
available  for  use  in  the  inversion  algorithm.  In  addition,  because  the  phase  is 
numerically  constructed  from  its  derivative,  the  spikes  in  the  latter  must  be 
selectively  removed  using  a  local  mean  amplitude  comparison  in  order  for  the 
algorithm  to  be  stable.  We  also  note  that  the  zero  crossings  in  Figure  68  occur  at 
different  depths  for  different  frequencies;  therefore,  for  this  profile,  the  phase 
monotinicity  is  a  function  of  frequency.  The  averaged  (over  depth  and  frequency) 
curves  for  the  amplitude  and  phase  derivative  are  illustrated  in  Figures  69  and  70, 
respectively.  Note  that  it  is  these  averaged  functions  that  will  be  used  in  an 
inversion. 

Earlier  in  the  discussion,  we  postulated  that  a  reconstruction  could  be 
attempted  over  a  compact  region  in  depth  if  the  phase  were  locally  monotonic. 
However,  we  first  need  to  consider  how  band-limiting,  in  particular  the  absence  of 
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Figure  67. 


Figure  68. 


AMPLITUDE  vs.  DEPTH 


Amplitude  of  DWB  kernel  for  10  source  frequencies  generated  using 
the  NORDA  2A  profile.  The  source  depth  is  .5  km  and  the  range  is 
30  km. 
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PHASE  DERIVATIVE  vs.  DEPTH 


OEPTH  (km) 


Derivative  of  the  phase  with  respect  to  depth  for  10  source 
frequencies  generated  using  the  NORDA  2A  profile.  The  source 
depth  is  .5  km  and  the  range  is  30  km. 
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AVERAGE  AMPUTUOE  vs.  DEPTH 


o 


Figure  69.  Frequency-averaged  amplitude  of  DWB  kernel  generated  using  the 
NORDA  2A  profile.  The  source  depth  is  .5  km  and  the  range  is  30 
km. 


AVERAGE  PHASE  DERIVATIVE  vs.  DEPTH 


Figure  70.  Frequency  and  depth-averaged  derivative  of  the  phase  with  respect 
to  depth  generated  using  the  NORDA  2A  profile.  The  source  depth 
is  .5  km  and  the  range  is  30  km. 
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zero  frequency,  affects  the  recovery  of  a  given  perturbation.  Also  of  importance  in 
the  reconstruction  process  is  the  effect  of  using  numerically  constructed  averaged 
quantities.  For  these  reasons,  we  next  investigate  the  effects  of  phase  averaging 
and  band-limiting  on  the  Fourier  reconstruction  of  two  different  perturbations. 

We  first  look  at  our  ability  to  reconstruct  a  simple  triangular  perturbation 
(see  Figure  71,  the  true  perturbation  has  a  maximum  value  of  1  at  100  m).  The 
perturbation  is  centered  at  100  m  and  its  extent  is  200  meters.  The  entire  region 
over  which  the  inversion  is  performed  is  assumed  to  be  700  meters.  Because  the 
perturbation  is  a  real  positive  function,  we  know  from  Fourier  analysis  that  zero 
frequency  is  required  to  fully  reconstruct  the  amplitude  of  the  perturbation. 
However,  very  low  frequencies  are  unavailable  in  an  actual  inversion,  and  so  first 
we  illustrate  the  effect  that  systematically  removing  the  lowest  frequency  values 
has  on  the  perturbation  recovery.  If  the  sampling  interval  is  taken  to  be  .684  m 
(1024  samples),  then  the  sampling  interval  in  frequency  (where  c0  is  1500  m/s)  is 
~1.0s_1.  If  we  generate  the  Fourier  transform  of  the  triangular  function,  zero  the 
low  frequency  components  of  the  transform  one  sample  at  a  time,  and  then  inverse 
transform,  we  obtain  the  curves  in  Figure  71.  The  actual  function  has  a  maximum 
amplitude  of  1.  Successively  smaller  amplitudes  represent  the  reconstruction  with 
one  less  frequency  component.  As  the  near-zero  components  are  discarded  the 
maximum  amplitude  of  the  function  decreases  correspondingly  and  the  function  is 
shifted  as  a  result  of  the  removal  of  the  zero  frequency  component. 
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LOW-CUT  FILTERED  PERTURBATION 
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Figure  71.  Effect  of  low-cut  filtering  on  the  Fourier  reconstruction  of  a  positive 
triangular  perturbation. 


We  next  investigate  the  effect  of  numerically  constructing  the  phase  from  the 
phase  derivative.  The  first  step  is  to  obtain  the  transform  of  the  triangular 
perturbation.  Next,  we  calculate  the  derivative  of  the  phase  of  the  transform 
kernel  and  from  this  function  we  integrate  to  find  the  unwrapped  phase.  If  the 
phase  derivative  is  a  smooth  function  of  depth,  then  the  previous  step  is  stable. 
However,  as  discussed  above,  the  phase  derivative  (e.g.  Figure  68)  has  large  local 
spikes  and  must  be  smoothed  before  calculation  of  the  phase  and  this  smoothing 
process  has  a  significant  effect  on  the  Fourier  inversion. 

Illustrated  in  Figure  72,  are  the  perturbation  reconstructions  obtained  from 
averaging  over  0  to  250  points  in  50  m  increments.  It  is  important  to  note  that 
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Figure  72.  Effect  of  phase-averaging  on  the  Fourier  reconstruction  of  a  positive 
triangular  perturbation. 
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these  curves  axe  the  real  paxt  of  the  function  obtained  from  the  inverse  transform. 
The  averaging  process  has  destroyed  the  original  symmetry  of  the  transform;  the 
spectrum  no  longer  corresponds  to  that  of  a  real  function.  Although  in  the  actual 
inversion  process  we  are  not  averaging  the  phase  of  the  data  (as  done  here),  but 
instead,  averaging  the  Fourier  phase,  the  effect  is  the  same;  components  needed  to 
reconstruct  the  perturbation  are  in  error. 

We  now  investigate  the  effects  that  band-limiting  and  phase  averaging  have 
on  the  recovery  ~f  a  second  perturbation.  The  perturbation  has  the  same  extent 
and  amplitude  of  the  previous  example,  but  now  is  placed  evenly  about  zero.  The 
effect  of  band-limiting  is  illustrated  in  Figure  73.  Because  the  DC  component  is 
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LOW-CUT  FILTERED  PERTURBATION 


Figure  73.  Effect  of  low-cut  filtering  on  the  Fou.ier  reconstruction  of  a 
symmetric  (about  zero)  perturbation. 
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Figure  74.  Effect  of  phase-averaging  on  the  Fourier  reconstruction  of  a 
symmetric  perturbation. 
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identically  zero,  elimination  of  the  near-zero  components  (0-4  Hz)  results  in  only 
half  the  amplitude  loss  as  compared  to  Figure  71. 

The  results  obtained  from  averaging  the  phase  of  the  spectrum  are  illustrated 
in  Figure  74.  The  reconstructions  are  more  consistent  as  compared  to  Figure  72; 
the  structure  is  well-identified  for  any  average  filter  length. 

From  the  examples  discussed  above,  it  appears  that  the  second  type  of 
perturbation  would  be  a  good  choice  to  try  in  the  inversion  algorithm.  However, 
no  consistent  reconstructions  could  be  obtained  for  either  type  of  perturbation  in 
the  NORDA  2A  environment.  In  fact,  although  we  will  briefly  discuss  the  deep 
ocean  examples  (test  cases  6  and  7,  chapter  4),  no  reconstructions  were  obtained 
for  any  profiles  studied. 

Several  explanations  for  the  failure  of  the  algorithm  can  be  offered.  First  is 
the  behavior  of  the  amplitude  and  phase  with  frequency.  Because  these  functions 
vary  with  frequency,  only  averages  are  available  for  use  in  the  reconstruction.  As 
previously  discussed,  the  effect  of  averaging,  in  particular,  phase  averaging  has  a 
strong  effect  on  the  c  p  t  of  the  Fourier  transform.  Second,  the  phase  is  not 
monotonic;  the  transformation  from  depth  to  phase  is  not  one-to-one.  Therefore, 
there  exists  a  neighborhood  (of  unknown  extent)  about  the  stationary  phase  point 
in  which  the  inversion  formulation  is  invalid. 

Although  we  are  unable  to  recover  a  perturbation  using  this  inversion 
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technique,  we  nonetheless  investigate  the  behavior  of  the  Born  kernel  for  the  deep 
ocean  profiles  (examples  6  and  7,  chapter  4)  in  order  to  gain  additional  insight. 
Illustrated  in  Figures  75,  76  and  77  are  the  transmission  loss  of  the  field  and  the 
amplitude  and  phase  derivative  of  the  Bom  kernel,  respectively,  constructed  using 
the  parameters  of  test  case  6.  Because  this  is  a  strongly  refractive  profile,  the 
structure  in  the  three  figures  is  sharply  defined.  The  energy  is  almost  purely 
refractive.  The  amplitude  and  phase  derivative  curves  for  10  frequencies  (5  to  50 
Hz)  at  a  range  of  30  km  are  illustrated  in  Figures  78  and  79  with  their  respective 
averages  in  Figures  80  and  81.  As  compared  to  the  Figures  67  and  68,  these  curves 
are  nearly  identical,  except  at  the  localized  spikes  and  low  frequencies.  In  this 
case,  the  stationary  points  consistently  occur  within  a  small  region  at  a  depth  of  2 
km.  However,  at  a  range  of  30  km  and  a  receiver  depth  of  500  m,  we  are  in  a 
shadow  zone.  From  the  modeling  results  of  chapter  4,  we  recall  that  the  DWB  and 
DWR  solutions  were  unable  to  give  a  valid  approximation  to  the  true  solution. 
Unfortunately,  if  we  move  farther  out  in  range  the  amplitude  and  phase  derivative 
behavior  again  varies  with  wavenumber.  Therefore,  although  the  DWB  kernel  is 
stable  (in  the  sense  that  the  Fourier  inversion  is  valid)  over  the  upper  region  of  the 
waveguide,  the  lack  of  field  structure  in  this  region  makes  a  perturbation 
reconstruction  impossible. 

In  the  last  example,  we  investigate  the  DWB  kernel  when  the  waveguide  is 
strongly  refractive  and  well-illuminated;  test  case  7.  The  transmission  loss  and  the 
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Figure  75. 


Transmission  loss  calculated  using  the  profile  and  source-receiver 
geometry  of  test  case  G.  The  source  frequency  is  hO  Hz. 


Amplitude  of  the  DWB  kerne!  calculated  using  the  profile  and 
source-receiver  geometry  of  test  case  6.  Note  the  trapping  of  energy 
in  narrow  bands  due  to  the  source  placement. 
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Figure  77.  Derivative  of  the  phase  with  respect  to  depth  calculated  using  the 
profile  and  source-receiver  geometry  of  test  case  6.  As  compared  to 
Figure  66,  there  are  broad  regions  in  depth  for  which  the  phase  is 
monotonic. 
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Figure  78. 


Figure  79. 


AMPLITUDE  VS.  depth 


0.0  2.0  «.o 
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Amplitude  of  DWB  kernel  for  10  source  frequencies  generated  using 
the  profile  and  source-receiver  geometry  of  test  case  6.  The  range  is 
30  km. 


PHASE  DERIVATIVE  vs.  DEPTH 


DEPTH  (km) 


Derivative  of  the  phase  with  respect  to  depth  for  10  source 
frequencies  generated  using  the  profile  and  source-receiver  geometry 
of  test  case  6.  The  range  is  30  km. 
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AVERAGE  AMPLITUDE  vs.  DEPTH 


Figure  80.  Frequency-averaged  amplitude  of  DWB  kernel  generated  using  the 
profile  and  source-receiver  geometry  of  test  case  6.  The  range  is  30 
km. 


AVERAGE  PHASE  DERIVATIVE  vs.  DEPTH 
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Figure  81.  Frequency  and  depth-averaged  derivative  of  the  phase  with  respect 
to  depth  generated  using  the  profile  and  source-receiver  geometry  of 

test  case  6.  The  range  is  30  km. 
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amplitude  and  phase  derivative  of  the  DWB  kernel  are  illustrated  in  Figures  82,  83 
and  84.  Although  the  structure  is  regular  as  compared  to  the  NORDA  2A  test 
case,  we  see  from  the  phase  derivative  that  there  are  many  zero  crossings  over  the 
depth  of  the  waveguide.  As  illustrated  by  the  amplitude  and  phase  derivative 
curves  versus  wavenumber,  (Figures  85  and  86)  and  their  averages  (Figures  87  and 
88),  the  behavior  of  the  DWB  kernel  varies  with  wavenumber.  In  Figure  88,  we 
observe  two  stationary  points  and,  in  addition,  a  higher  order  stationary  point  in 
the  vicinity  of  the  source.  Therefore,  for  this  source  geometry,  we  have  a  well- 
illuminated  waveguide,  but  the  points  of  stationarity  render  the  inversion  algorithm 
invalid. 


5.4  Summary  of  the  Inverse  Problem 

The  refractive  nature  of  propagation  in  the  ocean  waveguide  makes  a  wide¬ 
band  Fourier-like  inversion  for  a  single  source-receiver,  invalid.  In  regions  where 
the  waveguide  is  well-illuminated,  there  are  multiple  turning  points  (where  the 
Jacobian  goes  through  a  zero)  and  the  transformation  relating  the  scattered  field 
data  and  a  depth-dependent  perturbation  is  no  longer  1  to  1.  If  the  energy  is 
trapped  within  a  narrow  beam  due  to  source  placement,  then  the  regions  in  which 
the  phase  is  monotonic  coincide  with  shadow  zones  in  the  waveguide;  in  these 
regions,  the  DWB  and  DWR  approximations  yield  inaccurate  results. 


Figure  82.  Transmission  loss  calculated  using  the  profile  and  source-receiver 
geometry  of  test  case  7.  The  source  frequency  is  50  Hz. 
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Figure  83.  Amplitude  of  the  DWB  kernel  calculated  using  the  profile  and 
source-receiver  geometry  of  test  case  7.  Note  the  trapping  of  energy 
near  the  channel  axis. 


(km) 


Figure  84.  Derivative  of  the  phase  with  respect  to  depth  calculated  using  the 
profile  and  source-receiver  geometry  of  test  case  7.  As  in  Figure  66, 
there  are  many  zero  crossings.  In  this  case  the  turning  points  are 
due  to  refraction  within  the  waveguide. 
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0.0  2.0  4.0 


DEPTH  (km) 

Figure  85.  Amplitude  of  DWB  kernel  for  10  source  frequencies  generated  using 
the  profile  and  source-receiver  geometry  of  test  case  7.  The  range  is 
30  km. 
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Figure  86.  Derivative  of  the  phase  with  respect  to  depth  for  10  source 
frequencies  generated  usinj  the  profile  and  source-receiver  geometry 
of  test  case  7.  The  range  is  30  km. 
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AMPLITUDE  vs.  DEPTH 


Figure  87.  Frequency-averaged  amplitude  of  DWB  kernel  generated  using  the 
profile  and  source-receiver  geometry  of  test  case  7.  The  range  is  30 
km. 


PHASE  DERIVATIVE  vs.  DEPTH 


Figure  88.  Frequency  and  depth-averaged  derivative  of  the  phase  with  respect 
to  depth  generated  using  the  profile  and  source-receiver  geometry  of 
test  case  7.  The  range  is  30  km. 
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As  an  alternative  to  a  wavenumber  inversion,  we  could  formulate  the  inverse 
given  multiple  receiver  configurations,  either  horizontal  or  vertical  arrays,  for  a 
single  frequency  source.  For  example,  a  closely-spaced  horizontal  array  would 
permit  the  tracking  of  stationary  points  in  the  phase  of  the  kernel,  (which  are 
related  to  the  eigenray  arrivals)  and  therefore,  differentiation  of  the  arrivals. 
Recovery  of  a  profile  perturbation  could  be  attempted  by  first  evaluating  the  depth 
integral  at  all  receiver  ranges  using  the  method  of  stationary  phase.  If  we  then 
assume  that  the  magnitude  of  the  perturbation  at  the  stationary  point  for  a  given 
depth  is  the  same  at  all  receivers,  the  perturbation  at  each  depth  can  be  recovered 
using  a  matrix  inversion. 
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CONCLUSIONS 


The  Born  and  Rytov  approximations  to  the  transmitted  field,  in  which  the 
direction  of  propagation  is  in  range  and  the  perturbation  is  in  depth,  behave  in 
much  the  same  way  as  the  solutions  in  a  layered  earth.  In  a  constant  velocity 
earth,  the  depth-dependent  wave  functions  are  have  the  same  functional  form  for 
the  true  and  background  fields;  therefore,  the  errors  are  simple  to  derive  and 
interpret.  In  fact,  the  Born  and  Rytov  mode  coefficients  can  be  projected  from  the 
solutions  to  give  the  exact  error  in  the  approximations  for  a  single  mode.  As 
verified  by  the  constant  velocity  example,  the  validity  of  the  Born  solution  depends 
on  the  size  and  extent  of  the  perturbation;  the  scattered  field  must  be  less  than  the 
incident  field  in  order  to  obtain  an  accurate  approximation.  The  Rytov  solution, 
however,  in  insensitive  to  the  distance  traveled  in  the  perturbed  medium.  The 
accuracy  of  the  approximation  instead  depends  on  the  gradient  of  the  ratio  of  the 
scattered  to  incident  fields. 

The  behavior  of  the  solutions  in  a  layered  waveguide  is  a  bit  more  difficult  to 
quantize.  The  true  and  background  eigenfunctions  no  longer  have  the  same  form 
and  the  projected  mode  coefficients  are  difficult  to  interpret.  However,  through 
numerical  examples,  we  found  the  mode  functions  can  differ  by  up  to  7%;  at  this 
limit,  accurate  representations  of  the  true  field  will  be  obtained  for  ranges  less  than 
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25  km. 

The  information  obtained  from  the  two-layer  test  cases  is  used  to  estimate  the 
magnitude  and  extent  of  perturbations  permissible  when  the  background  sound 
speed  depends  on  depth.  The  accuracy  of  the  approximations  depends  both  on  the 
profile  and  the  source-receiver  configuration.  If  the  receiver  is  in  a  shadow  zone, 
the  correct  field  structure  is  not  well  approximated.  In  general,  however,  if  the 
receiver  is  well-  illuminated,  the  behavior  of  the  DWB  and  DWR  fields  (for  a 
depth-dependent  perturbation)  follows  that  predicted  in  the  two  layer  case. 
Accordingly,  the  extent  to  which  the  background  profile  can  deviate  from  the  true 
profile  in  an  inhomogeneous  environment  can  be  estimated;  therefore  indicating  the 
size  and  extent  of  profile  perturbatioi  yhich  can  be  recovered  in  an  inverse 
problem  using  either  the  DWB  and  DWR  approximations  to  the  scattered  field. 

The  inverse  algorithm  was  formulated  as  a  Fourier  inversion  assuming  the 
phase  75*  the  DWB  kernel  and  the  wavenumber  to  be  Fourier  transform  variables. 
Because  the  kernel  is  written  as  a  modal  sum,  no  explicit  expressions  for  the 
amplitude  and  phase  were  available;  therefore,  the  analysis  was  completely 
numerical.  Although  the  gross  wavenumber  dependence  was  determined  and 
numerically  removed  from  the  amplitude  and  phase,  the  functions  were  necessarily 
averaged  over  wavenumber  to  smooth  any  local  dependence.  When  the  averaging 
process  is  applied  to  the  phase,  the  information  content  of  the  Fourier  constituents 
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is  distorted;  in  effect  the  phase  function  is  no  longer  contins  the  correct  information 
on  the  perturbation.  In  addition,  the  phase  was  not  monotonic  in  regions  in  which 
the  waveguide  was  well-illuminated  and  so  in  these  regions  the  inverse  algorithm 
was  no  longer  valid.  As  a  result  of  these  limitations  (plus  band-limiting),  no 
perturbations  were  recovered  using  either  the  DWB  or  DWR  scattered  field. 
Therefore,  although  the  inversion  could  be  posed  using  olher  parameter 
configurations,  we  conclude  that  for  a  single  source-receiver  geometry,  a  wide-band 
Fourier  inversion  is  not  a  valid  formulation  for  recovering  depth-dependent 
perturbation  in  the  deep-ocean  waveguide. 
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APPENDIX  A: 

Derivation  of  the  Normal  Mode  Representation  of  the  Field 

The  integral  representation  of  the  sound  field  of  a  point  source  contains  the 
normal  modes  (discrete  spectrum)  as  well  as  the  continuous  spectrum.  The 
discrete  spectrum  constitutes  the  main  contribution  since  the  normal  modes  are 
only  weakly  attenuated,  whereas  the  continuous  spectrum  rapidly  attenuates  with 
range.  If  we  are  in  the  far  field  of  the  source,  then,  to  a  good  approximation,  only 
the  normal  modes  contribute  significantly  to  the  pressure  field.  In  this  appendix, 
the  modal  field  representation  is  derived  for  a  specific  set  of  boundary  conditions. 

We  begin  with  the  inhomogeneous  Helmholtz  equation 

V2  +  kln20(z)  P 0(k0,r, r0,z,z0)  =  (A1) 

in  cylindrical  coordinates  for  an  arbitrary  source  point,  ( r0,z0 ).  The  boundary 
conditions  are  chosen  to  be 

Q 

Po(k0,r,r0,z,z0 )  |  z=0  =  0,  ~Q^Po{koiriroizizo)  I  z=h  —  0 

where  h  is  the  width  of  the  ocean  waveguide. 

The  homogeneous  form  of  Eq.  (A.l)  is  separable  and  the  pressure  field  of  a 
wave  outgoing  from  a  source  located  at  r0  can  be  written  as 
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P0{ko, r,r0, z,z0)  =  C(z0,()H{2)((r)J0((r0)i/>(z,{)  r  >  r0  (A.2.a) 

Po{koJJo,z,zo)  =  C{z0,t)H(y[Zr0)Jo(tr)i)(z,t)  r  <  r0  (A.2.b) 

The  time  dependence  was  chosen  to  be  exp(iwf).  In  the  first  equation  H\?](£r) 

constitutes  a  right-traveling  wave.  The  second  equation  contains  both  right  and 
left-traveling  components,  but  at  large  r  the  right-going  wave  is  only  seen  in  its 
effect  due  to  the  coupling  with  Eq.  (A.2.a).  The  range-dependent  solution,  which  I 
will  call  G(£r  |  £r0)  for  simplicity,  is  derived  in  Appendix  B.  The  depth  dependent 
solution,  i>{z,£),  satisfies 

-77  +  (kWo(z)  ~  £2)  =  0  (A. 3) 

dz 

with  boundary  conditions 

=0  =  °- 
2=0  z=h 

This  system  is  solvable  only  for  specific  values  of  the  separation  parameter,  £. 

Chosing  two  linearly  independent  solutions  of  Eq.  (A. 3)  to  be  rpi(z,£)  and 
rp2(z,  0  then  the  depth  dependent  solution  is  written 

0(*>  £)  =  C'lV'i  (*,£)  +  C2ip2[z,$),  (A. 4) 

where  Cy  and  C2  are  constants.  Eq.  (A. 4)  and  the  boundary  conditions  will  give 
the  relationship  between  the  coefficients  for  all 
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jMgjn) 

C'~  1  MO, (n) 

and  also  the  characteristic  equation  for  the  eigenvalues,  £n 


(A.5) 


The  solution  of  the  inhomogeneous  equation  can  now  be  represented  as  the  sum  of 


normal  modes 


Po(K,r,'o,z,z„)  =  E  C(2o,(n)S((„r  j 


(A.6) 


The  excitation  coefficients,  C(z0,£n),  are  found  by  substituting  Eq.  (A.6)  into  the 
inhomogeneous  equation,  (A.l). 

n 

+  [cf(*o»€«)Gf(€»r  I  +  k20nl{z)C{z0,Zn)G(tnr  \  (nr0)tf(z,(„) 


ti(r~ro)S(z-Zo) 


(A.7) 


Using  the  equation  for  G(£nr  |  fr»ro)> 


d2  .  1  d 


77  +  ~rTr  +  **  I  ^  =  ~6<<r-ro)  ’ 

dr  r  ar  tt  r 


Eq.  (A.7)  reduces  to 


E  C(zo,tn)’l’(z,tn)  =  T6(Z~Zo) 


(A.8) 


(A.9) 
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Multiplying  by  rp(z,$m),  integrating  over  the  depth  of  the  waveguide  and  using 


orthogonality  of  the  normal  modes  gives 


c(z„  e„)  =  t 


4  h 


f  ll>(z,(nMz,tn)  & 

"  r\ 


(A.10) 


Normalizing  the  ip's  in  such  a  way  that  the  integral  is  equal  to  1,  C(z0,£n)  is  equal 

to  —  ip{z0,£n)>  and  the  fmal  solution  is  written  as 
4 

Po(koiriroiz>zo)  ~  'y  ]  G(£nr  |  £nr0)ip(z,(in)ip{z0,£n) .  (A. 11) 


Substituting  the  G(£nr  \  £nr0)  from  Appendix  B,  the  incident  field  modal  solution 


P0{k0,r,r0,z,z0)  =  jYj  ^2)  r>r0 


Po{koiriroizizo)  —  4  iS  [£nro)Jo{£nr)xt}{zy£n)‘tt}{zoi£n)  r<-ro •  (A. 12) 
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APPENDIX  B: 

Derivation  of  the  Range  Dependent  Solution 


Bessel’s  equation  for  an  arbitrary  source  point,  r  =  r0,  is 


7T  +  G^r  I  £«r°)  =  • 

dr*  r  ar  7rr 


(B.l) 


The  zero  order  forward  Hankel  transform  is  defined  to  be 


OO 

gfa)  =  J  9{r)J0{Kr)rdr 


and  the  corresponding  inverse  transform  is 


r°° 

<7(0  =  J  g(kf)J0{krr)krdkr 


where  kr  is  the  radial  wavenumber.  Applying  a  Hankel  transform  to  Eq.  (B.l), 
integrating  by  parts  twice  and  using  the  relationships  between  the  zeroth  and  first 
order  Bessel  functions  results  in  the  transformed  Green’s  function 


~  _  2 i  d0{kTr0) 

~  ”  k2-  ' 


(B.2) 


Inverse  transforming  gives 


G{Znr  |6»r„)  =  -  — / 


2 i  f°°  Jo{Kro)Jo{Kr) 


k2r  ~  Zn2 


kr  dkT 


(B.3) 


Using  the  relationship  between  Bessel  and  Hankel  functions 
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Jo{krro)Mkrr)  =  7 [M1)(^r'')  +  M2)(V)1  (B.4) 

4 


to  rewrite  the  integrand  and  the  relation 

BSH-Kr)  =  -HVHKr)  (B.5) 

to  reduce  the  four  integrals  to  two,  gives 


0° 

G{Znr  |  £nro)  =  ““  f 


H%\krr0)H$\krr)krdkr 


2tt  J  k2r  -  £„2 


— OO 

00 


-if 


—00 


^{kTr)H^{krr0)kTdkr 

k2r  -  e„2 


(B.6) 


If  we  now  change  the  integration  along  the  real  line  into  a  contour  integral  by 
making  kr  complex,  the  evaluation  of  these  integrals  is  accomplished  using  the 
residue  theorem.  If  the  time  dependence  is  chosen  to  be  e,u,i,  then  the  form  of  an 
outgoing  wave  is  e  **'r.  Using  the  asymptotic  form  of  the  Hankel  functions,  these 
integrals  can  be  rewritten  as 


G(£nr  |  (nr0) - 2~(rro)  1/2  /  7 - 7Tdkr 

*  C  kr  ~ 

-  v(rr°]~m  l^JiK 


(B-7) 


The  asymptotic  forms  are  used  simply  to  make  the  pole  evaluation  more 
transparent.  In  the  first  integral  of  Eq.  (B.7),  r  +  r0  is  always  greater  than  0; 
therefore,  if  kT  is  positive  imaginary,  we  have  convergence  in  the  upper  half  plane. 
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We  have  two  poles  at  which  to  evaluate  the  integrand.  Evaluation  at  the  positive 
pole  results  in  an  incoming  wave  and  so  this  solution  is  discarded.  Evaluation  at 
the  negative  pole  gives  the  outgoing  wave  contribution.  Asymptotically,  the  first 
integral  becomes 


-4-(rr„)-1/2e'i{'|r+r-), 

n 


(B.8) 


and  therefore,  the  solution  is 

(B.9) 

The  phase  of  the  second  integral  looks  like  kr(r  —  r0 )  and  must  be  evaluated 
depending  on  whether  r  is  greater  than  or  less  than  r0.  For  r>r0,  we  have 
convergence  in  the  upper  half  plane.  Evaluation  at  kT  =  —  £„  gives  the  outgoing 
wave  and  the  solution  is 


±W{i„r0)H<-V(inr).  (B.10) 

For  r<r0,  we  have  convergence  in  the  lower  half  plane.  Evaluation  at  the  positive 
pole  gives  the  outgoing  wave  solution 

Summing  Eqs.  (B.9)  and  (B.10)  and  using  Eq.  (B.4),  we  obtain  the  solution  for  r 


greater  r0. 
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GUnrl^r.)  *  Y^(fi»r0)J5?l2)(£llr)  r>r0.  (B.12) 

Summing  Eqs.  (B.9)  and  (B.ll)  and  using  Eq.  (B.4),  we  obtain  the  solution  for  r 
less  than  r0 

C(f„r|enr„)  =  r<r..  (B.13) 

Eqs.  (B.12)  and  (B.13)  are  the  range-dependent  Green’s  function  for  a  point  source 
located  at  r  =  r0. 

Note  that  Eqs.  (B.12)  and  (B.13)  are  a  general  result  for  r0  =*=  0.  In  the 
special  (and  usual)  case  where  the  source  is  located  at  r0  =  0,  we  see  from  Eq. 
(B.12)  that  we  recover  the  usual  radial  Green’s  function 

<5(«„r|?„r»)  =  r.=0. 


(B.14) 
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APPENDIX  C: 

Derivation  of  the  Full  DWR  Solution 


The  DWR  approximation  is  derived  by  first  assuming  the  total  field  solution 
can  be  written  as 


P[Ks,ro,z,z0)  =  Po{K,r,ro>z,Zo)^k°'r'T°'Z'Z°]  , 


(C.l) 


where  P0{k0,r, r0,z,z0)  is  a  known  incident  field  and  $(k0,r,r0,z,z0)  is  a  complex 
phase  function.  Substituting  Eq.  (C.l)  into  the  Helmholtz'  equation,  (l),  and 
performing  the  differentiation  gives 


V2  +  |V2$  +  (V$)2  +  Jfc»»2(r,z)j  P0  +  2 VP0-V$ 


6{r-r0)6(z-z0) 


e"*.  (C.2) 


From  Chapter  2,  we  know  the  incident  field  satisfies 


V2  +  k2n2{z)^P0  =  - 


6[r-ro)6{z~zo) 


(C.3) 


Therefore,  by  evaluating  e  ®  at  the  source  ($  is  zero  at  the  source  because  the 
scattered  field  is  zero),  Eq.  (C.2)  can  be  reduced  to 


P0  V2$  +  2  WVV$  =  -P0  (V$)2  -  k2(n2(r,z)  -  n20(z))P0  .  (C.4) 


We  would  like  to  write  Eq.  (C.4)  in  the  form  of  a  Helmholtz  equation  and  solve  for 
$  using  our  Green’s  function,  P0.  Away  from  the  source,  we  can  set 
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P0  V2$  +  2  VP0-V$  = 


V2  +  k\nl(z ) 


(Po*) 


(C.5) 


so  using  Eq.  (C.5)  our  differential  equation,  (C.4),  becomes 


V2  +  k2n2(z) 


(P0*)  =  -P0  (V$)2  -  k\(n2(r,z)  -  nl(z))P0  .  (C.6) 


Therefore,  the  complex  phase  function,  $,  is  given  by  the  non-linear  integral 
equation 


$  = 


2irlci 


oo  h 


~f  f  Po(ko,r,r'>z,z')(n2(r,z)  -  n2(z))P0(k0,r' ,r0,z’ ,z0)r' dr' dz' 


+  ~~J  J />o(*0,r,r',z,z')[V$(r',z')]2 P0(k0,r’ ,rotz’ ,z0)r' dr' dz' .  (C.7) 


o  o 

oo  h 


0  0 


To  obtain  the  first  DWR  approximation,  we  expand  the  complex  phase  in 


powers  of  e 


$(ko,r,r0,z,z0)  =  <f>\  +  <£2  +  •  •  •  (C.8) 

and  use  Eq.  (4)  (the  expansion  on  the  difference  of  the  squares  of  the  index  of 
refraction).  We  substitute  the  expansions  into  Eq.  (C.7),  and  equate  powers  of  e. 

The  second  line  in  Eq.  (C.7)  is  assumed  to  be  order  e2  due  to  the  square  of  the 
phase  gradient.  This  suggests  that  the  field  must  vary  smoothly  in  order  for  the 
first  DWR  approximation  to  be  valid.  The  first  integral  has  one  term  of  order  e; 
therefore,  to  0(e),  the  complex  phase  is  approximated  as 
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${k»>riro7z>zo)  ~  Mko,r,r0,z,z0) 


'/ f  f>o(ko,r,r',z,z')nl(r',z')P0(k0,r',r0,z',z0)r'dr'dz'  (C.9) 


and  the  first  DWR  approximation  is 


d®1  =  p.t' 


(C.10) 


as  derived  from  the  first  DWB  in  Chapter  2. 
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APPENDIX  D: 

Evaluation  of  the  Range-Dependent  Solution 

As  defined  in  chapter  2,  the  range-dependent  solution  for  a  point  source 
located  at  r  =  0  is  given  by 

r 

I[n,m)  =  ffi2)($nr)f  {Zmr')r’ dr' 

O 

+  U(nr)f  ^(i„r')^Umr')r'dr- .  <DI) 

r 

For  clarity  in  the  evaluation,  the  first  integral  is  defined  to  be  Ix  and  the  second  to 
be  /2.  The  integrals  are  evaluated  differently  depending  on  whether  £„  =  or 

£n*£m- 

For  the  integrals  are  evaluated  as  follows.  We  know  J0{£nr)  and 

satisfy  Bessel’s  equation.  Choosing  x(r)  and  y(r)  to  represent  any 
Bessel  functions,  we  write 

x"  +  —X  +  i2  x  =  0 
r 

and 

y"  +  — y'  +  6n2  y  =  o. 

r 

Multiplying  the  first  equation  by  r  y,  the  second  by  r  x,  subtracting  the  two 
equations  and  combining  terms  gives 
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0  0  0  0 
yx  -  xy 


+  x'  y  -  y'  x  =  (£„2  -  £mz)rx  y 


or 


r(xy'-yx')]  =  £„2  -  £m2 


r  xy 


To  evaluate  the  first  integral  in  Eq.  (D.l),  we  set  J0(£nr)  =  x>  =  y  and 

integrate  from  0  to  r  to  obtain 

r 

Wm-aS.OflP’K.O  -  r' 


(D.2) 


m 


The  solution  of  the  second  integral  is  obtained  in  a  similar  manner. 
Multiplying  the  terms  by  their  respective  constants,  /x  and  I2  become 


'  -  tj 

MM 


* 2  -  2  i  2 


(n  -  e. 


Ur'H[2H(mrVMnr')  ~  £.rV1(f,r’)Sl!l(£„r') 

-  «„r'M21  (€»r ')«<*> ') 


(D.3) 


Summing  the  contributions  from  the  limit  at  r  gives 


=  /, 


+  Ii 


inf 


c  2  _  t  2 


(D.4) 


Substituting  the  relations 
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•M£„r)  =  {  [flf’fe.r)  +  flj?'(£„r) 
A(«,r)  =  |  [fll"  («„r)  +  M2)  (£.0 


(D.S) 


into  Eq.  (D.4)  and  recognizing  the  Wronskian 


1 
2  L 


fl?)(£.r)MI)(£.r)  -  M2'(£„r)M1)(£nr) 


2 1 


ntn?  ’ 


(D.6) 


we  obtain  the  result 


7(n,m) 


2t  HVHtmr) 

*  tn*  -  £m2  ' 


(D.7) 


Evaluation  of  Ix  at  zero  givc3  0.  Evaluation  of  /2  at  r  =  oo  is  zero  assuming  that 
in  and  have  a  small  negative  imaginary  part.  The  solution  then  for  in=*=im  is 
given  by  Eq.  (D.7). 

The  evaluation  for  in  =  im  is  somewhat  more  complicated.  Setting  im  —  in 
explicitly  in  the  integrals  gives 


r 

I{n,n)  =  H^][inr)  J  J„{inr')H^]  {inr  ')r'dr ' 

O 

OO 

+  f  fl?>(£.r-)flp>  (£.r>'*'. 

r 


(D.S) 


The  first  step  is  to  make  a  change  of  variables.  Setting  y=inr'  and  dr'  =  dy/i 
the  integrals  become 
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I  =  ^ ’,nr-f  J0[y)H{o){y)ydy  +  —  J  Hi2]{y)Hi2){y)ydy  .  (D.9) 

Zn  o  („r 

Again,  we  define  the  first  integral  as  Ix  and  the  second  as  /2  and  integrate  by 
parts.  To  evaluate  Ilt  we  write 


u  =  J0(y)H{2)(y)  du  =  -  J0{y)H^]  (y)  +  J1(y)Hi02)(y) 


dy 


dv  =  ydy 
to  obtain 


v  = 


h  = 


M2)  (e„r) 


in 


■^-4(s)M2)(»)  {"  -/  4rP.(sOM21  [«)  +  J\(y)li,':](y),  dy 


<d.io) 


Evaluating  the  first  term  at  £„r,  (the  expression  is  zero  at  y  =  0),  we  obtain 


|jjp,(e.0 


(D.ll) 


To  evaluate  the  remaining  integral,  /u,  in  Eq.  (D.10),  we  substitute  Eqs.  (D.5)  for 
J0(y )  and  ^i(y)  and  combine  terms.  The  result  is 


/n  =  - 


■f  [2^(y)M2)(y)  +  M1}(y)M2)(y)  +  M1’ (y)^,2)(y)]  dy.  (D-12) 
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Defining  the  first  term  as  Jin  and  using  Bessel’s  Equation 

y2M2)(y)  =  -y^;(2)(y)  -  y2#;,(2)(y) 

we  obtain 


(m  = 


h?KU')  r("r 

2$n2 

M2)(fnr) 


/  [yfl*.<2)(»)  +  y2J^;'(2)(sr)]  H-M(y)  iy 


2in  dV 


f-’jL 
J  dv 


jiL 


W(2)(y))5 


dy. 


Integrating  and  evaluating  at  the  limits,  / m  becomes 


'ill 


=  [M2)  (€«r) 


Writing  the  last  two  terms  as  /112,  using  the  relations 

rfa'Hy)  =  My)  +  <K(y) .  fl»2)M  =  My)  -  <Y0(y) 
h[2] (y)  =  A(y)  -  >YM  ,  Mll(y)  =  Mv)  +  >T,(y) 


and  combining  terms  gives 


A12  - 


fll2'  (€»r) 


2f 


n  o 


f  [y2^(y)^o(y)  +  y2  ^(y)  ^(y) 


dy 


Using  Bessel’s  equation,  /112  becomes 


fll2)(£„r) 


yi(e.r)i  +  (v1U.r) 


The  second  integral  in  Eq.  (D.9)  is  found  using  the  same  techniques. 


(D.13) 


(D.14) 


(D.15) 


(D-16) 


(D.17) 


(D.18) 

but  the 
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evaluation  is  simplier  because  the  integrand  is  the  product  of  identical  Bessel 
functions.  The  solution,  found  by  integration  by  parts  and  Bessel’s  equation  is 


!■>  =  - 


+  (M2)  (c„>-)): 


(D.19) 


The  terms  contributing  to  the  solution,  I  =  Ix  +  /2,  are  in  Eqs.  (D.ll),  (D.15), 
(D.18)  and  (D.19).  Canceling  like  terms,  we  obtain 


1=  - 


MM 


L  J  4 


.  *!?>(«.••) r-  r  12 

+  - - -  J \  (£nr)j  + 


2  (£„r)  r2  r  12 

+  -—7--—  Y^M 

4  L 


(D.20) 


Writing  tf<„2>({„r)  =  J.( f„r)  -  iY. (€„r)  in  the  second  term  and  recombining 


terms  gives 


I =  - 


*  .  ^4,!,(f»r)r2 


[M"  («.r)  (^‘'(e.r))’]  (D.21) 


"I  * 

where  *  denotes  complex  conjugate.  Now,  H[1^  (£„r)  =  Hi2'1  (£„  r),  so  for  real 


eigenvalues 


/  =  - 


Si‘>((nr) 


—  [M2)  (e~0 


2 


(D-22) 


or  recognizing  the  Wronskian 


-  M1)(e»r)M2) (enr)  =  -■ 


(D.23) 


we  obtain  the  solution  when  as 
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if  M2)(gnO 

7T  £n 


(D  .24) 


Therefore,  evaluation  of  the  range  integral  for  a  depth-dependent  profile  results  in 
the  solutions 


2 

j.  •>  +  1  ‘  ^ 

*■  (J  -  (»2 

ir  M2)  (f-r)  ,  , 


(D.25.a) 


(D.25.b) 


It  is  interesting  to  note  that  Eq.  (D.25.b)  also  follows  by  applying  l’Hospital’s  rule 
to  Eq.  (D.25.a),  in  which  case  the  result  is  not  confined  to  real  eigenvalues. 
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APPENDIX  E  : 

Modal  Decomposition  of  the  Perturbative  Solutions 

In  this  Appendix,  we  discuss  the  modal  decomposition  of  the  DWB  and  DWR 
solutions.  The  DWB  and  DWR  mode  coefficients  axe  first  derived  for  the  full  field 
solutions.  These  expressions  are  later  reduced  to  that  of  a  single  mode  calculation 
and  the  mode  coefficients  are  given  explicitly  for  the  test  cases  in  Chapter  3. 

The  DWB  solution  for  a  strictly  depth-dependent  profile  is  given  in  Eqs.  (18) 
and  the  DWR  solution  in  Eqs.  (19).  Multiplying  both  sides  of  the  DWB  solution 
by  the  eigenfunctions,  t/>(z,£/),  and  integrating  over  the  width  of  the  waveguide,  we 
obtain 

h 

f  r,  0, z, z0 )  r{j(z,  &)  dz 

O 

h 

=  f  [Po(h,r,0,z,zo)  +  Pan(ko,r,0,z,zo)  ip(z,£i)  dz  £„  =  (E-l  a) 
0 

and 

h 

f  P[B}  (ko,r,0,z,zo)t/>(z,$i)  dz 

0 

h 

=  /  [P0(ko,r,0,z,zo)  +  P)nrn(tco,r,0,z,zo)  r/>(z,$,)  dz  (E'Lb) 

O 

Performing  the  same  operations  on  the  DWR  solution  gives 
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n 

f  P[R)  {^o,r,0,z,zo)  xp{z,£i)  dz 


i 

=  / 


f  ~  \ 


P0{K^O,z,z0)ex  p 


rl>(z,ii)dz  £n  =  £m  (E.2.a) 


and 


n 

f  P{R)  (h,r,0,z,zo)  ip(z,$i)  dz 


t 

-5 


P0{.Kir,Q,z,Zo)ex  P 


snm 


(E.2.b) 


Because  the  eigenfunctions  are  orthonormal  over  the  width  of  the  waveguide,  the 
dz  integrals  are  1  if  /  =  n  and  zero  otherwise.  Eqs.  (E.l)  and  (E.2)  are  the  general 
expressions  used  in  determining  the  DWB  and  DWR  mode  coefficients.  For 
simplicity,  the  expressions  will  now  be  evaluated  for  a  single  mode  for  the  test 
cases  discussed  in  Chapter  3. 

The  first  example  is  of  a  constant  velocity  waveguide  in  which  the 
perturbation  spans  the  entire  waveguide.  Eqs.  (E.l.b)  and  (E.2.b)  do  not 
contribute  because  the  perturbation  is  constant.  In  addition,  we  are  finding 
solutions  for  only  1  mode.  From  Eq.  (26),  the  eigenfunctions  are  sines;  replacing 
the  ij)'s  with  sines  (see  also  Eq.  (18. b))  in  Eq.  (E.l. a)  and  combining  terms,  the 
single  mode  Bom  expression  becomes 
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1/2  * 


\  f  Pn(B){ko,r,0,z,zo)sinX1zdz  =  j#i2)(£ir) 


2]  7  w7r2°'i  f  • 

-  sm (— )J  sm(-)  - 


.  ,TTZ.  ,r^nl  .  /7T2n  .  /7T2,  , 

;m(  — — ) - - — sin  —  smf —  <22 

*2A'  2^  v2/i'  ^  2h ' 


(E.3) 


where  Aj  has  been  replaced  by  — —  .  Integrating  the  right  hand  side  yields 

2  h 


l1/2  /•* 

T  /  pii,Z’l(*:o.r.0.;!.2o)sin(ff) 


-7^1 


zVA;^  n ! 

2^i 


(E.4) 


Dividing  through  by  the  range  solution,  “#i^(£ir)>  the  resulting  Born  mode 

4 

coefficient,  Cu j,  in  a  constant  velocity  waveguide  is  given  by 


-  i  sto<if > 1  +  it 


i>^"i 


(E.5) 


Because  the  Rytov  solution  is  essentially  an  exponentiated  Born,  to  avoid  the 
simple  analysis,  we  will  simply  state  the  result.  The  Rytov  mode  coefficient,  CiR 


Ctp  = 


.  *z0.  [  irk20nx 

Sm  ~2/T  eX^  ~W 


The  second  example  for  which  the  first  mode  coefficient  is  derived  is  a  two-layer 


waveguide.  The  velocity  is  constant  in  each  layer  and  the  background  profile  a 
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constant  over  the  entire  waveguide.  Again,  Eqs.  (E.l.b)  and  (E.2.b)  are  not 
evaluated  because  we  are  interested  in  single  mode  analytic  solutions. 

In  this  case,  the  mode  functions  are  not  the  same  for  the  true  and  background 
profiles.  The  true  mode  functions  axe  given  in  Eq.  (38)  and  the  background  modes 
are  sines.  The  Bom  and  Rytov  approximations,  although  constructed  from  the 
guess  modes,  are  multiplied  by  the  true  mode  functions  to  obtain  the  mode 
coefficients. 

From  Eqs.  (E.l.a),  (38)  and  (40. a),  the  Born  projection  equation  is 


h 

f  P[B)  {k0,r,Q,z,z0)Tl)x {z,£)dz  =  -^•sin(^-)M2)(^i>')v/^7 


- - —f  sin(-^-)  sinAu  dz  +  f  sin(-^-)  cosA0l(/i-z)  dz 

smXnZi  •  2n  *  2fx 


ik?0rnx 


7T2„ 


+  — 


•  (*Zl\ 

sm<lT) 


7T 


v^T 


cosAo^/i-zx)  f  1  .  ,nz^  .  %  ,  ,  f  .  ,ttzs  ,  ,L  s 

- — -  I  sm(— ) smAu  dz  +  f  sin(— )cosA 0l[h-z)  dz 

sinAuZj  J  2h  2  h 


(E.7) 


Integrating,  combining  terms  and  dividing  by  the  range  solution  yields 

iklrnx 


cib  =  v^T  sin(^-) 


i  + 


il 

h 


•  (nZl\ 

sm(— j 


7T 


(E.8) 
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where  y/~A~i  is  given  by 


\f~A~i  =  Ai  16h2sin(-^-)  cosAo^fe-Zi) 


r  \  7TZi 

2hXn2XoltanXol(h — z\)  2hX0i  XncotXnZi-\-Tr  A0i  An  cot(  ) 


2h 


jr2  -  4h2Aol2  7r2  -  4A2An2 


Accordingly,  the  Rytov  is 


(E.9) 


(E.IO) 


